K Zhang

D Sciacca

A Coati

R Bernard

Y Borensztein

P Diener

B Grandidier

I Lefebvre

M Derivaz

C Pirri

G Prévot

Resolving the structure of the striped Ge layer on Ag(111): Ag 2 Ge surface alloy with alternate fcc and hcp domains

Two-dimensional (2D) honeycomb lattices beyond graphene, such as germanene, promise new physical properties such as quantum spin Hall effect. While there have been many claims of growth of germanene, the lack of precise structural characterization of the epitaxial layers synthesized hinders further research. The striped layer formed by Ge deposition on Ag(111) has been recently ascribed as a stretched germanene layer. Using surface X-ray diffraction and density functional theory calculations, we demonstrate that it corresponds in fact to a Ag2Ge surface alloy with an atomic density 6.45% higher than the Ag(111) atomic density. The overall structure is formed by stripes associated with a face-centered cubic top-layer alignment, alternating with stripes associated with an hexagonal-close-packed top-layer alignment, in great analogy with the (22 × √3) Au(111) reconstruction.

Introduction

After the first experimental observation of the growth of Ge on Ag(111) in 1999 [START_REF] Oughaddou | Ge/Ag(111): SURFACE ALLOY OF A SEMICONDUCTOR ON A METAL[END_REF][START_REF] Oughaddou | Ge/Ag(111) Semiconductoron-Metal Growth: Formation of an Ag2Ge Surface Alloy[END_REF], the germanium/silver interface has been subject to a renewed interest since the discovery of germanene, the graphene-like allotrope for Ge atoms [START_REF] Cahangirov | Two-and One-Dimensional Honeycomb Structures of Silicon and Germanium[END_REF][START_REF] Derivaz | Continuous Germanene Layer on Al(111)[END_REF]. Ge structures formed upon deposition at room temperature on Ag(111) have been initially observed by low energy electron diffraction (LEED) and Auger electron spectroscopy (AES) [START_REF] Oughaddou | Ge/Ag(111): SURFACE ALLOY OF A SEMICONDUCTOR ON A METAL[END_REF]. Two ordered superstructures were found, corresponding to different Ge coverages. The first structure, corresponding to a �√3 × √3�𝑅𝑅30° reconstruction related to the Ag(111) surface (hereafter named as �√3 × √3�), was assigned to a Ge coverage Ge θ =1/3 ML, with respect to the atomic density of a Ag(111) plane. From core-level photoemission spectroscopy (XPS), a model of Ag2Ge surface alloy, where every three Ag atom is replaced with a Ge atom was proposed for the �√3 × √3� structure and confirmed from angle resolved photoemission spectroscopy (ARPES) [START_REF] Golias | Surface Electronic Bands of Submonolayer Ge on Ag(111)[END_REF]. The second structure, corresponding to a (7 × 7) reconstruction, was assumed to correspond to a complete Ge monolayer. This reconstruction was further recognized as a c�√3 × 7� reconstruction and a model of Ge tetramers on top of a Ag(111) plane was proposed for its structure, based on scanning tunneling microscopy (STM) [START_REF] Oughaddou | Germanium Adsorption on Ag(111): An AES-LEED and STM Study[END_REF]. In contradiction with this model, from real-time STM observations of the growth, we have recently proposed that this latter structure is also a surface alloy with a coverage of 0.6 ± 0.1 ML [START_REF] Zhang | Growth of Germanium-Silver Surface Alloys Followed by in Situ Scanning Tunneling Microscopy: Absence of Germanene Formation[END_REF].

The �√3 × √3�𝑅𝑅30° structure obtained after deposition at room temperature was found to evolve upon annealing at 473 K [START_REF] Wang | Broken Symmetry Induced Band Splitting in the Ag 2 Ge Surface Alloy on Ag(111)[END_REF]. A striped pattern, with ridges and valleys, appears in STM images, and LEED diffraction spots are split into satellite spots. Thus, this structure was shown to deviate from the simple Ag2Ge surface alloy and assigned to a �6√3 × √3�𝑅𝑅30° reconstruction. This structure can also be obtained by evaporating Ge at 600K [START_REF] Liu | Role of Atomic Interaction in Electronic Hybridization in Two-Dimensional Ag 2 Ge Nanosheets[END_REF]. The striped pattern has been attributed to the relaxation of the compression force induced by the difference in the radii between Ge atoms and Ag atoms after Ge insertion [START_REF] Liu | Role of Atomic Interaction in Electronic Hybridization in Two-Dimensional Ag 2 Ge Nanosheets[END_REF].

The Ge/Ag(111) system has been recently revisited and contradictory observations have been reported. On the one hand, the striped phase has been interpreted as a highly stretched pure germanene layer, with 23% and 12% stretch in the direction parallel and perpendicular to the stripes [START_REF] Lin | Single-Layer Dual Germanene Phases on Ag(111)[END_REF]. On the other hand, in a recent study, the LEED diagram corresponding to satellite spots around diffraction conditions of a �√3 × √3� reconstruction, and previously attributed to the striped phase, has been interpreted as a �19√3 × 19√3�𝑅𝑅30° reconstruction, corresponding to a Ag2Ge surface alloy contracted by 5% with respect to the Ag(111) surface [START_REF] Chiniwar | Substrate-Mediated Umklapp Scattering at the Incommensurate Interface of a Monatomic Alloy Layer[END_REF].

In any event, in spite of several studies devoted to Ge/Ag(111), the precise atomic structure of the reconstruction formed at low coverage and for deposition above room temperature or after annealing is unknown. In particular, there is no consensus about the formation of a layer of germanene or of a Ag-Ge alloy. Surface X-ray diffraction (SXRD) is very well adapted for elucidating the structure of ordered surface reconstructions, and we have recently used it to successfully determine the precise atomic positions for the various silicene epitaxial layers on Ag(111) [START_REF] Curcella | Determining the Atomic Structure of the ( 4 × 4 ) Silicene Layer on Ag(111) by Combined Grazing-Incidence x-Ray Diffraction Measurements and First-Principles Calculations[END_REF][START_REF] Curcella | Structure and Stability of Silicene on Ag(111) Reconstructions from Grazing Incidence x-Ray Diffraction and Density Functional Theory[END_REF].

In this paper, we present a SXRD study of the �√3 × √3� Ge/Ag reconstruction. We show that the ordered structure associated with the satellite spots around �√3 × √3� diffraction conditions can be precisely indexed as a c(31 × √3) reconstruction, and thus corresponds to the striped phase observed in STM images. Comparison with DFT calculations allows us to determine the precise atomic structure of the surface, which corresponds to a Ag2Ge surface alloy with an atomic density 6.45% higher than the Ag(111) atomic density. The overall structure is formed by stripes associated with a face-centered cubic (fcc) top-layer alignment, alternating with stripes associated with an hexagonal-close-packed (hcp) top-layer alignment, in great analogy with the (22 × √3) Au(111) reconstruction.

Methods

SXRD experiments were performed at the SIXS beamline of SOLEIL synchrotron. The Ag(111) sample was prepared by repeated cycles of Ar+ sputtering and annealing at T=750 K. Ge was evaporated in the diffraction chamber from a crucible using a Knudsen cell with a sample kept at ∼ 420 K. The flux was estimated to 0.5 ML/h. The sample was analyzed with 18.46 keV X-rays at a grazing incidence angle of 0.2°. Scattered X-rays were detected with a X-ray Pixel chips with Adaptive Dynamics hybrid pixel detector [START_REF] Dawiec | Real-Time Control of the Beam Attenuation with XPAD Hybrid Pixel Detector[END_REF]. Diffracted intensity was measured by performing a series of rocking scans around diffraction. We used the "BINoculars" software to produce threedimensional (3D) intensity data in the reciprocal space from the raw data [START_REF] Roobol | BINoculars : Data Reduction and Analysis Software for Two-Dimensional Detectors in Surface X-Ray Diffraction[END_REF]. The intensity was further integrated along the direction parallel to the surface to obtain the structure factors. For this purpose, the data were fitted with the product of a lorentzian lineshape in one direction with a gaussian lineshape convolved with a door lineshape in the other direction, using a home-made software, and the fitted function has been analytically integrated. We finally obtained a set of 2493 Theoretical computations of the minimal energy configuration were done within the density functional theory (DFT). Calculations for the core electrons were performed using the projector augmented wave (PAW) method [START_REF] Blöchl | Projector Augmented-Wave Method[END_REF][START_REF] Kresse | From Ultrasoft Pseudopotentials to the Projector Augmented-Wave Method[END_REF], as implemented in the VASP code [START_REF] Kresse | Efficiency of Ab-Initio Total Energy Calculations for Metals and Semiconductors Using a Plane-Wave Basis Set[END_REF][START_REF] Kresse | Efficient Iterative Schemes for Ab Initio Total-Energy Calculations Using a Plane-Wave Basis Set[END_REF]. The Perdew-Burke-Ernzerhof functional [START_REF] Perdew | Generalized Gradient Approximation Made Simple[END_REF] was used to describe the exchange-correlation functional. The plane wave basis set was restricted to a cut-off of 300 eV. The Ag(111) substrate was modeled by a (31 × √3) nine-layer slab. Similar to the procedure adopted in ref. [START_REF] Golias | Surface Electronic Bands of Submonolayer Ge on Ag(111)[END_REF][START_REF] Curcella | Structure and Stability of Silicene on Ag(111) Reconstructions from Grazing Incidence x-Ray Diffraction and Density Functional Theory[END_REF], the bottom two layers were kept fixed at the equilibrium theoretical positions. The 10 th layer on top is modeled as a Ag2Ge surface alloy following the SXRD results. The vacuum region is 9 Å thick. The entire system, Ag2Ge alloy on silver, was fully relaxed by a conjugate gradient method until the forces acting on each atom were less than 0.01 eV/Å. The convergence criterion for self-consistent field calculations in energy was chosen of 10 -4 eV. A 1 × 8 × 1 k-point mesh was used to sample the Brillouin zone. Figure 1a shows the diffracted intensity for in-plane conditions (𝑙𝑙 = 0.12) after evaporation of Ge θ ≈1/3 ML at T = 420 K. Satellite spots appear around diffraction spots corresponding to the crystal truncation rods (CTR) of the substrate (integer values of ℎ and 𝑘𝑘 indices). Other sets of spots appear around fractional values of (ℎ, 𝑘𝑘), i.e. diffraction conditions corresponding to a

Experimental results

Fig

(√3 × √3) reconstruction (for example, ℎ = 1 3 , 𝑘𝑘 = 1 3
). A detailed view of the diffracted intensity around ( 13 , 1 3 , 0.12) and (1, 0, 0.12) is shown in Fig. 1b and1c.

q (Å -1 )
These spots are aligned along three directions, indicated by white dotted lines in Fig. 1b and1c.

These <110> directions are equivalent due to the p3m1 symmetry of the substrate surface. Thus, the reconstruction should correspond to the three possible orientations of a similar reconstruction.

We can exclude the fact that all spots belong to a single domain of trigonal or hexagonal symmetry since it would imply a large number of extinctions which are not observed in the present experiments.

A profile along the A-A' line shown in Fig. 1c is drawn in Fig. 2. The central spot corresponds to the CTR and 3 orders of diffraction from each side are visible for the satellite spots. The spacing between the spots is ∆𝑞𝑞 = 0.1412 ± 0.0002 Å -1 , indicating a periodicity of 44.5 Å, i.e, 15.4 times the Ag-Ag interatomic distance 𝑎𝑎 𝐴𝐴𝐴𝐴 = 2.889 Å. Thus, as expected from previous works [START_REF] Wang | Broken Symmetry Induced Band Splitting in the Ag 2 Ge Surface Alloy on Ag(111)[END_REF], the measured phase appears as a periodic modulation of a (√3 × √3) reconstruction, along the <110> directions. A precise analysis of the diffraction diagram show that the unit cell is a rectangular centered unit cell. As 15.4 is close to 31/2, it can be described as a rectangular c(31 × √3)

reconstruction. The value of 15.4𝑎𝑎 𝐴𝐴𝐴𝐴 is in good agreement with the periodicity of the striped pattern previously observed in STM images [START_REF] Lin | Single-Layer Dual Germanene Phases on Ag(111)[END_REF]. Thus, the diffraction pattern measured by SXRD corresponds to the striped pattern previously observed on STM images. On the contrary, we can exclude that this structure corresponds to a �19√3 × 19√3�𝑅𝑅30° reconstruction as proposed in ref. [START_REF] Chiniwar | Substrate-Mediated Umklapp Scattering at the Incommensurate Interface of a Monatomic Alloy Layer[END_REF].

The intensity and the shape of the spots in Fig. 1b and 1c are modulated by the elongated X-ray beam footprint on the sample, in the real space, and by the intrinsic width of the diffraction pattern, in the reciprocal space. For spots around CTRs, the intrinsic width is mainly given by the finite size of the striped domains and by the fluctuation of the striped phase periodicity. The first contribution is the same for all the spots whereas the second one increases linearly with the distance of the satellite spot to the substrate spot. The width of the satellite spots is measured to ∆𝑞𝑞 = 1.2 10 -3 + 6 10 -4 𝑛𝑛 Å -1 where n is the diffraction order. This shows that the periodicity is very well defined, and the lateral size of the domains is of the order of 0.5 µm. (see Fig. 1). This observation is in good agreement with a chemical ordering between Ag and Ge Figure 4 shows the configuration for the alloyed structure obtained after relaxation. The 𝑦𝑦(𝑥𝑥)

DFT computations

profile indicates that the atomic positions of surface atoms oscillate between fcc (at 𝑦𝑦 = 1/3) and hcp sites (at 𝑦𝑦 = 1/6). The 𝑧𝑧(𝑥𝑥) profiles in Fig. 4c shows that on average, Ge atoms are located around 0.1Å below the Ag atoms. The 𝑧𝑧(𝑥𝑥) profiles for Ge and Ag atoms display a periodic oscillation with a double frequency as compared with 𝑦𝑦(𝑥𝑥). Starting from 0, the odd and even minima for 𝑧𝑧 correspond to atoms in fcc and hcp positions respectively, whereas maxima correspond to atoms in bridge position. This gives rise to the periodic striped pattern observed in STM images [START_REF] Wang | Broken Symmetry Induced Band Splitting in the Ag 2 Ge Surface Alloy on Ag(111)[END_REF]. The calculations indeed indicate an apparent periodicity of 22.4 Å with to 5.5 Ge atoms per stripe, in good agreement with the 6√3 periodicity found by STM [START_REF] Wang | Broken Symmetry Induced Band Splitting in the Ag 2 Ge Surface Alloy on Ag(111)[END_REF]. The value of the buckling is 0.12Å, whereas the buckling measured by STM is approximately 0.2 Å [START_REF] Wang | Broken Symmetry Induced Band Splitting in the Ag 2 Ge Surface Alloy on Ag(111)[END_REF]. The Ge-Ag interatomic distances are slightly smaller for Ge atoms in bridge position (2.662 Å) than for Ge atoms in fcc position (2.722 Å). The value of the relaxed atomic positions is given in the Supplemental Material [START_REF]for Fitting Parameters and Comparison between Experimental and Simulated Structure Factors[END_REF]. This undulation of the atomic positions between the fcc sites and the hcp sites is strongly analogue to the (22 × √3) Au(111) reconstruction, for which the surface atomic density is 23/22 times higher than the one of a Au(111) bulk plane, and which displays also a similar striped structure [START_REF] Narasimhan | Elastic Stress Domains and the Herringbone Reconstruction on Au(111)[END_REF].

Comparison with SXRD and discussion

From this atomic configuration relaxed by DFT, we have computed the simulated structure factors where 𝑁𝑁 pts = 2493 is the number of experimental structure factors, 𝑁𝑁 par is the number of free parameters and 𝜎𝜎 exp is the experimental uncertainty, which takes into account the statistical uncertainty given by the number of counted photons and an overall 10% uncertainty estimated from the comparison of symmetryrelated structure factors.

Fig. 5 shows the comparison between the experimental and simulated structure factors 𝐹𝐹 for inplane measurements and along selected rods. Since 𝐹𝐹(𝐻𝐻, 𝐾𝐾, 0) = 𝐹𝐹(-𝐻𝐻, -𝐾𝐾, 0) and 𝐹𝐹(𝐻𝐻, 𝐾𝐾, 0) = 𝐹𝐹(-𝐻𝐻, 𝐾𝐾, 0) due to the mirror symmetry of the unit cell, only the quadrant corresponding to 𝐻𝐻 > 0, 𝐾𝐾 > 0 has been drawn. The complete comparison of 62 non-equivalent reconstruction rods is given in Fig. S1 and the Debye Waller parameters used for the fit are given in Table S1 [START_REF]for Fitting Parameters and Comparison between Experimental and Simulated Structure Factors[END_REF]. There is an outstanding agreement between experiments and simulations, corresponding to a small value 𝜒𝜒 2 = 1.90. As can be seen, the simulation reproduces very well both satellite rods near integer values of (ℎ, 𝑘𝑘) (Fig. 5b-e) and satellite rods near fractional values of (ℎ, 𝑘𝑘) (Fig. 5f-i). The first ones display intense variations at specific integer values of 𝐿𝐿, for example, near 𝐿𝐿 = 1 and near 𝐿𝐿 = 4 in Fig. 5b and5e, near 𝐿𝐿 = 2 in Fig. 5c, or near 𝐿𝐿 = 0 and near 𝐿𝐿 = 3 in Fig. 5d. These positions are close to the Bragg diffraction conditions of the Ag(111) crystal. These intense variations can be attributed to the periodic elastic relaxations that penetrate in the bulk [START_REF] Croset | Elastic Relaxations and Interactions for Vicinal and Self-Organized Surfaces: Role of Crystalline Anisotropy[END_REF][START_REF] Prévot | Direct Observation of Elastic Displacement Modes by Grazing-Incidence X-Ray Diffraction[END_REF].

Such relaxations are induced by the 6.45% misfit between the outermost surface layer and the substrate. Thus, the fact that they are nicely reproduced indicates that DFT simulations give a very accurate description of the interaction between the surface layer and the substrate, even if the number of layers free to relax has been limited to 7 in the DFT calculations.

In order to estimate the uncertainties related to the atomic positions, we have let all atomic positions free to move in the first two planes, starting from the values obtained by DFT. This leads to a huge number of free parameters (𝑁𝑁 par =111) and to a reduction of 𝜒𝜒 2 down to 0.9. The atomic positions are not strongly modified. Displacements from the initial positions are less than 0.2 Å for Ag and Ge atoms in the first layer and less than 0.1 Å for atoms in the second layer, with rms displacements of 0.08 Å, 0.08 Å and 0.05 Å in the first layer along the x, y, and z directions respectively. Similar values are 0.03 Å, 0.03 Å, 0.06 Å for atoms in the second plane.

For the sake of completeness, we have also compared the experimental structure factors to the simulated structure factors obtained within a model of a germanene epitaxial layer. For that purpose, a model with 22 Ge hexagons along the long side of the mesh has been used, so that the Ge atomic density is twice the one of the alloyed model. The atomic positions obtained after relaxation by DFT are drawn in Fig. S2. As previously, we have set identical DW factors for all Ge atoms, for all Ag atoms of the second plane and for all other Ag atoms. We find a strong deviation between experiments and simulations as 𝜒𝜒 2 = 26.5. Therefore, the germanene layer cannot account the experimental observation and this model must clearly be rejected.

We have also computed the adsorption energy for Ge atoms in the two models:

𝐸𝐸 ad = (𝐸𝐸 Ge-Ag -𝐸𝐸 Ag -∆𝑁𝑁 Ag 𝐸𝐸 Ag bulk )/𝑁𝑁 Ge Where 𝐸𝐸 Ge-Ag is the total energy of the considered model, 𝐸𝐸 Ag is the energy of a slab of same lateral size without Ge, 𝐸𝐸 Ag bulk = -2.72 eV is the bulk cohesive energy of Ag atoms, 𝑁𝑁 Ge the number of Ge atoms in the slab, and ∆𝑁𝑁 Ag is the difference of number of Ag atoms between the considered Ge/Ag model and the bare Ag slab.

𝐸𝐸 ad is lower for the Ag2Ge surface alloy (-4.56 eV) than for the germanene layer (-4.41 eV).

Taking as a reference the Ge bulk cohesive energy (-4.45 eV), the germanene layer is less stable than the Ge bulk whereas the Ag2Ge surface alloy is more stable. This indicates that, from a thermodynamic point of view, wetting of the Ag surface with a Ag2Ge alloy is favored, whereas whereas a germanene layer would be unstable and should dewet to form 3D Ge crystallistes. Thus, there is a better thermodynamic stability for the Ag2Ge surface alloy than for the germanene layer.

The atomic structure of the striped phase presents a great analogy with the (22 × √3) Au(111)

reconstruction. This reconstruction corresponds to a surface atomic density 23/22 times higher than the one of a Au(111) bulk plane [START_REF] Narasimhan | Elastic Stress Domains and the Herringbone Reconstruction on Au(111)[END_REF]. For the Au surface atoms which are less coordinated than the bulk ones, the strength of the bonds is higher than in the bulk. Their equilibrium interatomic distance would thus be lower than in the bulk. This results in a large tensile stress for surface atoms. For Au(111), the system relaxes in order to reduce the interatomic distances for surface atoms, i.e., by increasing the surface atomic density. The energy gain due to the decay of the interatomic distances at the surface exceeds the energy cost of the regions where atoms occupy bridge positions instead of three-fold coordinated hcp or fcc positions. Unlike Au(111), Ag(111)

does not spontaneously reconstruct. However, when a Ge atom replaces an Ag atom of the surface, this increases the absolute value of the tensile surface stress since the atomic radius of Ge (1.25 Å) is less than the one of Ag (1.60 Å). Above a critical coverage, the system relaxes in a configuration where the interatomic distances are decreased, at the cost of the creation of Shockley partial dislocations, i.e., discommensuration lines separating regions where atoms occupy fcc positions, from regions where atoms occupy hcp positions. Other systems have been reported to show a similar behavior such as Cu, Ag, Au films on Ru(0001) [START_REF] Pötschke | Interface Structure and Misfit Dislocations in Thin Cu Films on Ru[END_REF][START_REF] Ling | Herringbone and Triangular Patterns of Dislocations in Ag, Au, and AgAu Alloy Films on Ru[END_REF].

Conclusion

Using GIXD, we have unambiguously determined the structure of the striped pattern observed upon Ge deposition on Ag(111) at 420 K, which was controversially attributed to a pure germanene layer or to a Ge-Ag alloy. From a careful analysis of the Patterson map, leading to an initial structural model of the system, we obtained the model of the Ag2Ge surface alloy with an atomic density 6.45% higher than the one of Ag(111). The computed theoretical structure factors obtained from this model, relaxed by DFT, show an outstanding agreement with the experimental ones, demonstrating that this striped phase is actually a Ag2Ge surface alloy. The observed stripes correspond to Shockley partial dislocations that separate alternating fcc and hcp domains, presenting a strong analogy with the Au (22 × √3) reconstruction. A model of germanene can be ruled out since it does not fit the experiments and since it is thermodynamically less stable.

It has to be noted that such a determination was possible thanks to the very large set of acquired structure factors, which can now easily be done within a relatively short time thanks to 2D detectors. This shows how powerful DFT associated with GIXD can be for the analysis of complex surface structures. In particular, since DFT is shown here to give a very accurate description of the atomic relaxations, the comparison between GIXD and DFT is a criterium of choice for determining the validity of an atomic model of a surface structure.

  structure factors along 62 inequivalent reconstruction rods. The (ℎ, 𝑘𝑘, 𝑙𝑙) indices used for indexing a reflection in reciprocal space refer to the Ag(111)) surface basis (𝑎𝑎 0 = 𝑏𝑏 0 = 2.889 Å, 𝑐𝑐 0 = 7.075 Å, 𝛼𝛼 0 = 𝛽𝛽 0 = 90°, 𝛾𝛾 0 = 120°). The (𝐻𝐻, 𝐾𝐾, 𝐿𝐿) indices refer to the Ag(111) c(31 × √3) reconstruction basis (𝑎𝑎 = 89.545 Å, 𝑏𝑏 = 5.003 Å, 𝑐𝑐 = 7.075 Å, 𝛼𝛼 = 𝛽𝛽 = 𝛾𝛾 = 90°).

. 1 . 1 3 , 𝑘𝑘 = 1 3 )

 113 Fig. 1. Diffracted intensity for in-plane conditions (𝑙𝑙 = 0.12) after evaporation of Ge θ ≈1/3 ML.

Fig. 2 .

 2 Fig. 2. Variation of the diffracted intensity for in-plane conditions around the (ℎ = 1, 𝑘𝑘 = 0, 𝑙𝑙 =

Fig. 3 .( 1 × 1 )

 311 Fig. 3. Experimental Patterson map of the c(31 × √3) structure. As it is centered, only half of the unit cell has been drawn along 𝑥𝑥. (size of the Patterson map is 44.78 × 5.00 Å 2 ).

Fig. 4 .

 4 Fig. 4. Top view (a) and profiles (b, c) along 𝑥𝑥 for surface atoms (blue: Ge, grey: Ag). The scales

Fig. 5 .

 5 Fig. 5. (a) Comparison between experimental (red half-disks) and simulated (black half-disks) in-

  and compared them to the experimental ones. Only two scale factors (one for the in-plane set and one for the rods) and Debye-Waller (DW) factors along the three directions were used as free parameters. For simplicity, we have set identical DW factors for all Ge atoms, for all Ag surface atoms, for all Ag atoms of the second plane and for all other Ag atoms. Thus only 14 free parameters were used to fit the data. The agreement between experimental and simulated (𝐹𝐹 th ) structure factors is estimated by the value of 𝜒𝜒 2 = 1 𝑁𝑁 pts -𝑁𝑁 par ∑ � 𝐹𝐹 th -𝐹𝐹 exp 𝜎𝜎 exp � 𝑁𝑁 Pts 2

Acknowledgments

This study is financially supported the French National Research Agency (Germanene project ANR-17-CE09-0021-03). K.Z. is supported by the Chinese Scholarship Council (CSC contract 201808070070).

atoms. Indeed, if one third of the surface atoms are Ge atoms, the intensity associated with (√3 × √3) satellites should roughly scale as �𝑍𝑍 Ge + 2𝑍𝑍 Ag cos � From these observations, we propose that the striped phase has a density 33/31 times higher than the Ag(111) atomic density, and that the atomic positions undulate between the fcc sites and the hcp sites. Such structure would thus display a great analogy with the (22 × √3) Au(111) reconstruction for which the surface atomic density is 23/22 times higher than the one of a Au(111) bulk plane [START_REF] Narasimhan | Elastic Stress Domains and the Herringbone Reconstruction on Au(111)[END_REF]. Consequently, starting from this Ag2Ge alloy surface model, we have relaxed


Juan-Pablo Cerapio

Marion Perrier

Camille-Charlotte Balança

Pauline Gravelle

Fréderic Pont

Christel Devaud

Don-Marc Franchini

Virginie Féliu

Marie Tosolini

Carine Valle

Fréderic Lopez

Anne Quillet-Mary

Loic Ysebaert

Alejandra Martinez

Jean Pierre Delord

Maha Ayyoub

Camille Laurent

Jean-Jacques Fournie
email: jean-jacques.fournie@inserm.fr

Juan-Pablo Orcid

Cerapio

Phased differentiation of γδ T and T CD8 tumor-infiltrating lymphocytes revealed by single-cell transcriptomics of human cancers

Keywords: Gamma-delta lymphocyte, scRNAseq, differentiation

γδ T lymphocytes diverge from conventional T CD8 lymphocytes for ontogeny, homing, and antigen specificity, but whether their differentiation in tumors also deviates was unknown. Using innovative analyses of our original and ~150 published single-cell RNA sequencing datasets validated by phenotyping of human tumors and murine models, here we present the first high-resolution view of human γδ T cell differentiation in cancer. While γδ T lymphocytes prominently encompass TCRVγ9 cells more differentiated than T CD8 in healthy donor's blood, a different scenario is unveiled in tumors. Solid tumors and lymphomas are infiltrated by a majority of TCRVγnon9 γδ T cells which are quantitatively correlated and remarkably aligned with T CD8 for differentiation, exhaustion, gene expression profile, and response to immune checkpoint therapy. This cancer-wide association is critical for developing cancer immunotherapies.

Introduction

The presence and functional differentiation of cytotoxic tumor-infiltrating lymphocytes (TILs) is critical for the clinical outcome of cancer patients. Although T CD8 TILs are associated with favorable prognosis in all types of cancers, little is known about other cytotoxic TILs of γδ T lineage. Human γδ T cells are CD4 -CD8 -T lymphocytes expressing either TCRVδ2 or TCRVδ1/3 chains prominently paired with TCRVγ9 or TCRVγ2/3/4/5/8 (TCRVγnon9), respectively. On the one hand like T CD8 cells, circulating γδ T lymphocytes differentiate into naïve (Tn), central memory (Tcm), effector memory (Tem), and terminally differentiated effector memory CD45Ra+ (Temra) cells. [START_REF] Dieli | Differentiation of effector/memory Vdelta2 T cells and migratory routes in lymph nodes or inflammatory sites[END_REF] Similarly, both can infiltrate tumors and mediate cytotoxic functions important for antitumor immunity. On the other hand, γδ T lymphocytes differ from conventional T CD8 by tissue distribution, relative abundance, and a more pleiotropic immunobiology. [START_REF] Hayday | The innate biologies of adaptive antigen receptors[END_REF] Unlike T CD8, γδ T cells react to TCR-independent and TCR-dependent stimuli that differ notably from major histocompatibility complex (MHC)-restricted peptides. Peripheral TCRVγ9 γδ T lymphocytes are specific for phosphoantigen-sensing butyrophilins, while tissular TCRVγnon9 γδ T cells recognize more diversified antigens.

Whether T CD8 and γδ T lymphocytes infiltrating tumors are somehow associated remains unclear, but it is important to develop and monitor γδ-targeted cancer immunotherapy. [START_REF] Sebestyen | Translating gammadelta (gammadelta) T cells and their receptors into cancer cell therapies[END_REF] By deconvoluting 10,000 tumors from 55 cancer types, we and others have previously found that the abundance of TCRVγ9 TILs is unrelated to that of TCRαβ TILs, varies considerably across individuals and cancers, and is associated with better outcomes in cancers. [START_REF] Tosolini | Assessment of tumorinfiltrating TCRVγ9Vδ2γδlymphocyte abundance by deconvolution of human cancers microarrays[END_REF][START_REF] Meraviglia | Distinctive features of tumor-infiltrating gammadelta T lymphocytes in human colorectal cancer[END_REF][START_REF] Gentles | The prognostic landscape of genes and infiltrating immune cells across human cancers[END_REF] Likewise, the presence, function, and potential of TCRVγnon9 cells which represent another prominent subset of γδ TILs are demonstrated in different human tumors. [START_REF] Wu | An innate-like Vδ1+γδ T cell compartment in the human breast is associated with remission in triple-negative breast cancer[END_REF][START_REF] Mikulak | NKp46expressing human gut-resident intraepithelial Vdelta1 T cell subpopulation exhibits high antitumor activity against colorectal cancer[END_REF][START_REF] Foord | Characterization of ascites-and tumor-infiltrating gammadelta T cells reveals distinct repertoires and a beneficial role in ovarian cancer[END_REF] So far, however, neither the general landscape of both subsets of cytolytic γδ TILs nor their putative link with the T CD8 had been characterized both at the single cell level and at a cancer-wide scale.

To address this issue, the characterization of single-cell transcriptomes from γδ TILs requires to identify exhaustively and selectively these cells in multiple scRNAseq datasets from different human tumors. Although theoretically only the γδ T lymphocytes express a cell surface γδTCR, other cell subsets do actually express some level of TCRγ or TCRδ mRNAs. Both TCRαβ T CD8 lymphocytes and NK cells also express either of the TRDC, TRGC1, or TRGC2 genes, imposing the use of a multigene signature rather than single genes to reliably identify γδ T lymphocytes only. [START_REF] Pizzolato | Single-cell RNA sequencing unveils the shared and the distinct cytotoxic hallmarks of human TCRVdelta1 and TCRVdelta2 gammadelta T lymphocytes[END_REF] Here, we incepted a 'score and gate' strategy that significantly improves the digital purification of γδ T lymphocytes from scRNAseq datasets, and developed a cell pseudotime-based visualization of differentiation trajectories. With the γδ T and the T CD8 lymphocytes extracted from ~150 human samples of PBMC and tumor biopsies, the pseudotimed differentiation trajectories of these lineages revealed their strikingly contrasting hallmarks in blood and cancer.

Results

Circulating γδ T are more differentiated than T CD8 lymphocytes in healthy individual's blood

From seven scRNAseq datasets of healthy donor's PBMC and γδ T lymphocytes cell-sorted from PBMC, [START_REF] Pizzolato | Single-cell RNA sequencing unveils the shared and the distinct cytotoxic hallmarks of human TCRVdelta1 and TCRVdelta2 gammadelta T lymphocytes[END_REF] (n = 6,122) γδ T lymphocytes were digitally extracted by 'score and gate' (Supplemental Data Figure S1) (Methods). Then, a differentiation trajectory of these extracted cells was computed, and the trajectory's principal components were projected onto the linearized pseudotime axis to visualize a pseudotimed trajectory (Methods). The pseudotimed differentiation trajectory of peripheral blood γδ T lymphocytes comprised 14 milestones (Figure 1a) in which TCRVγ9 were prominent over TCRVγnon9 cells, as expected (not shown). This trajectory was anchored by expression of Tn genes including the lymph node homing receptor CCR7 [START_REF] Girard | HEVs, lymphatics and homeostatic immune cell trafficking in lymph nodes[END_REF] and the TCF1 differentiation transcription factor-encoding gene TCF7, [START_REF] Zhou | Differentiation and persistence of memory CD8(+) T cells depend on T cell factor 1[END_REF] which peaked in the M1, M2 milestones. Further M3-M7 displayed intermediate levels of Tn and Tem gene expression, a defining hallmark of Tcm lymphocytes. [START_REF] Sallusto | Central memory and effector memory T cell subsets: function, generation, and maintenance[END_REF][START_REF] Restifo | Lineage relationship of effector and memory T cells[END_REF] The Tcm M6 and M7 upregulated both Tn (CD27, CD28, SELL, IL7R) and tissue migration genes (CD44, CXCR3, CXCR4). The Tem milestones M10-M11 included activated cells (IFNG, CCL5), and cytotoxic Tem cells (NKG7, PRF1, GZMA, GZMB, GZMM, FASLG) in M10-M13. The maximal expression of FCGR3A (CD16 marker), PRF1, and GZM genes defined the Temra cells, [START_REF] Angelini | FcgammaRIII discriminates between two subsets of Vgamma9Vdelta2 effector cells with different responses and activation pathways[END_REF] in M14. The exhaustion transcription factor TOX [START_REF] Khan | TOX transcriptionally and epigenetically programs CD8(+) T cell exhaustion[END_REF] and immune checkpoint genes LAG3, TIGIT, and PDCD1 were weakly expressed in control γδ T cells (Figure 1b). When projected onto the trajectory, the scores for naïve (Tn) and effector (Teff) gene signatures evolved smoothly and inversely, reflecting the transcriptional progression along aging from naïve to Temra cells (Figure 1c). Projecting on this trajectory, several published T CD8 differentiation signatures (Table S1) confirmed the differentiation stage of each milestone (Supplemental Data Figure S2), whose ontogeny was validated by the dendrogram of milestones (Figure 1d).

With 15 milestones, the pseudotimed trajectory of the (n = 3,680) PBMC-derived T CD8 lymphocytes approximated that of γδ T cells (Figure 1e), showing cells of stem central memory profile (IL7R, CD28) [START_REF] Gattinoni | T memory stem cells in health and disease[END_REF] in Tcm M4, and cells with progenitor exhausted-like (Tpex) profile (CD95, GZMK, TOX, TIGIT) [START_REF] Galletti | Two subsets of stem-like CD8(+) memory T cell progenitors with distinct fate commitments in humans[END_REF] in M6. The activated Tem cells (IFNG) upregulated PDCD1 and LAG3 in M11, while the T CD8 Temra (FCGR3A) cells peaked at final M15 instead of final M14 in PBMC (Figure 1f). These identifications were supported by projecting the above Tn, Teff, and other published signatures (Figure 1g, Supplemental Data Figure S3) and by the dendrogram of T CD8 milestones (Figure 1h).

Furthermore, (n = 137) γδ T and (n = 725) T CD8 lymphocytes were extracted from our CITEseq dataset of healthy donor's PBMC immunostained for differentiation markers, [START_REF] Pont | Single-cell virtual cytometer allows user-friendly and versatile analysis and visualization of multimodal single cell RNAseq datasets[END_REF] and were mapped onto the above trajectories. Their differentiation stages determined by the cell surface CD45RA and CD62L immunophenotypes were consistent with those assigned by gene signatures (Supplemental Data Figure S4). Reflecting distinct biologies and antigen specificities, however, the proportions of cells at each differentiation stage differed between the lineages and across individuals: γδ T were prominently Tem, while Tn predominated more often among T CD8 (Figures 1i, Supplemental Data Figure 4c). Therefore, in healthy adults PBMC, the γδ T are generally more advanced than T CD8 lymphocytes on their respective differentiation trajectories.

Convergent differentiations of γδ T and T CD8 lymphocytes infiltrating tumors

Then, to compare the differentiation hallmarks of γδ T and T CD8 lymphocytes infiltrating tumors (TILs), we performed 3ʹ chemistry scRNAseq of malignant lesions from (n = 3) splenic marginal zone lymphoma (SMZL) (GSE175785), and (n = 4) ovarian carcinoma (OVCAR) (GSE148162) cancer patients. Additional scRNAseq datasets from (n = 129) additional human tumors of various cancer types were also downloaded (summarized in Table S2), normalized, and integrated for comparison (Methods). Altogether, a total of (n = 4,648) γδ T and (n = 45,174) T CD8 TILs were extracted from a total of 136 tumors (Tables S3-S6), and analyzed as above (Supplemental data).

The tumors richest in γδ T TILs corresponded to the richest in T CD8 TILs (Fischer Exact p = .002), and the richest in TCRVγ9 were also the richest in TCRVγnon9 cells (Fischer Exact p = .0001). The trajectories of each lineage of TILs superimposed to those of PBMC (Figure 2a), but some milestones enriched in cancer revealed that both lineages of TILs (Figure 2b) were more differentiated than the PBMC (Figure 1i). The cell counts of γδ T and T CD8 TILs were strongly correlated (Pearson r = 0.7-0.99) in most cancer types except in lung carcinoma (LUAD) and melanoma (MEL). In Hodgkin's lymphoma (HL) and head and neck carcinoma (HNSCC), the viral status of patients unveiled a higher correlation of the γδ T and T CD8 TILs counts in the EBV-positive and HPV-positive, respectively, than in their negative counterparts (Supplemental Data Figure 5).

Classifying the tumor infiltration extents by TIL cell counts (Methods) indicated that the OVCAR and HNSCC were frequently inflamed (I) while the SMZL and hepatosplenic carcinoma (HCC) had classic (C) infiltrates, contrasting with the immune deserts (D) more frequent in many solid tumors (e.g. LUAD, MEL). The γδ T/T CD8 TIL ratio decreased from desert to inflamed infiltrates (mean D: 24%, C:10%, I: 5%) (Figure 2c, Table S2b). Identifying the Tex and Ttrm cells in all tumors (Methods) showed their respective rates were correlated between the two lineages (Pearson r = 0.8 in both cases) (Figure 2d). Confirming recent reports, [START_REF] Wu | An innate-like Vδ1+γδ T cell compartment in the human breast is associated with remission in triple-negative breast cancer[END_REF][START_REF] Mikulak | NKp46expressing human gut-resident intraepithelial Vdelta1 T cell subpopulation exhibits high antitumor activity against colorectal cancer[END_REF][START_REF] Foord | Characterization of ascites-and tumor-infiltrating gammadelta T cells reveals distinct repertoires and a beneficial role in ovarian cancer[END_REF] the majority of γδ TILs cancer-wide was of TCRVγnon-9 subtype, reflecting the extra-lymphoid tropism of this subset (Figure 2e, Table S2b). Combining trajectory and signature analyses revealed that, like Tex T CD8 TILs, all the Tex γδ TILs displayed Tem or Tcm 2e,f). Comparatively to their non-exhausted counterparts, these γδ T and T CD8 Tex cells shared a signature of 14 up-regulated genes involved in T cell exhaustion (Supplemental Data Figure 6a,b). These lymphocytes were also more glycolytic and more abundant among, but not restricted to, the Ttrm cells (Figure 2g, h, Supplemental Data Fig. 6 c). In the non-exhausted compartment, the few Tn and early Tcm γδ TILs presented gene signatures of resting unexpanded T CD8 TIL clones (tumor singletons), [START_REF] Wu | Peripheral T cell expansion predicts tumour infiltration and clinical response[END_REF] and the most mature Tem γδ T TILs had signatures of peripherally expanded large T CD8 clones 21 associated with response to immune checkpoint blockade (ICB) therapy [START_REF] Wu | Peripheral T cell expansion predicts tumour infiltration and clinical response[END_REF][START_REF] Yost | Clonal replacement of tumor-specific T cells following PD-1 blockade[END_REF] (Supplemental Data Figure 6c,d). As a result, the single-cell composition of all cancer samples showed that the differentiated states of γδ T and T CD8 were far more coherent in TILs than in PBMC.

Ex vivo immunophenotypes of γδ T and T CD8 TILs are correlated

To validate these conclusions, we obtained tumor biopsies from five patients with untreated ovarian cancer (OVCAR, n = 2), head and neck squamous cell carcinoma (HNSCC, n = 1) and cervical cancer (CC, n = 2), isolated their TILs, and immunostained the cell surface markers of γδ T and T CD8 lymphocytes for flow cytometry analysis.

Consistent with the scRNAseq data, both lineages of TILs from these lesions were prominently CCR7 -CD45RA -Tem cells (Figure 3a). In these epithelial cancers, the tumor antigenspecific Tex CD8 TILs associated with response to ICB are highly cytotoxic PD-1 hi TIM-3 + cells co-expressing CD39 and Ttrm markers. [START_REF] Balanca | Dual relief of T-lymphocyte proliferation and effector function underlies response to PD-1 blockade in epithelial malignancies[END_REF] So, we analyzed these features at the cell surface of the γδ T and T CD8 TILs from the above biopsies, and found sizable fractions of PD-1 hi TIM-3 + cells and of PD-1 hi -CD39 + cells among both lineages (Figure 3b). Despite the previously depicted variability of TILs across tumors (Figure 2b), the samples with high Tex rates in T CD8 TILs had also high Tex rates in γδ T TILs (Methods) and the same applied for samples with low Tex rates (Figure 3b). In all biopsies, both the CD8 T and γδ T lymphocytes encompassed PD-1 hi TIM-3 + and PD-1 hi CD39 + cells, and across these samples, their respective rates of T CD8 and γδ T TILs were correlated (Figure 3c).

We further analyzed (n = 12) mice from a syngeneic murine model of HNSCC obtained by orthotopically implanting the lung carcinoma TC1 cells retrovirally transduced to express HPV16 E6/E7 antigens. [START_REF] Balanca | Dual relief of T-lymphocyte proliferation and effector function underlies response to PD-1 blockade in epithelial malignancies[END_REF] The TILs were isolated from established tumors for immunostaining and flow cytometry analyses. These assays showed the emergence of Tex γδ T and T CD8 TILs at the same time points of growth development, and the rates of PD-1 hi TIM-3 + Tex cells in these γδ T and in T CD8 TILs were correlated (Figure 3d). In these mice, of the 16 possible co-expression patterns of the PD-1, TIGIT, CTLA-4, and TIM-3 immune checkpoints, [START_REF] Balanca | Dual relief of T-lymphocyte proliferation and effector function underlies response to PD-1 blockade in epithelial malignancies[END_REF] the same 5 phenotypic patterns predominated among both γδ T and T CD8 TILs (Figure 3e). Collectively, these ex vivo data confirmed the correlated differentiation patterns of γδ T and T CD8 TILs.

Both T CD8 and γδ T TILs are impacted by ICB therapy

The convergent differentiation of intratumoral γδ T and T CD8 lymphocytes reported above prompted us to investigate whether γδ T TILs are as affected as T CD8 TILs to immune checkpoint blocking (ICB) therapy. We downloaded two published scRNAseq datasets of 21 tumors from (n = 11) BCC [START_REF] Yost | Clonal replacement of tumor-specific T cells following PD-1 blockade[END_REF] and (n = 10) MEL [START_REF] Sade-Feldman | Defining T cell states associated with response to checkpoint immunotherapy in melanoma[END_REF] patients-derived tumor biopsies sampled before and after ICB treatment. Their γδ T and T CD8 TILs were extracted as above, mapped on their respective differentiation trajectory, and analyzed with regard to patient's clinical response, tumor infiltration pattern, and infiltrate's reactivity to treatment (post-treatment TIL counts increase) (Methods).

In all samples, the T CD8 outnumbered the γδ T cells (Figure 4a). The pre-treated BCC lesions were more infiltrated than those of MEL (mean: 295 versus 86 TILs, respectively), and were more reactive to treatment (BCC: 8/11 reactive, mean counts: 306 (pre)/785 (post), MEL: 4/10 reactive, mean counts: 73 (pre)/149 (post)) (Figure 4a). In most BCC, the treatment increased the cell counts for both lineages and differentiation stages, including Tex, without affecting their relative rates of Tn, Tcm, Tem, and Temra cells (Figure 4b,c). Although the prominently desert infiltrates of MEL hampered further conclusions, the BCC results showed that ICB therapy strengthens the association of tumor infiltrating T CD8 and γδ T lymphocytes (Pearson r = 0.77 pre-Tt, 0.95 post-Tt). Among these latter, however, although similar counts of TCRVγ9 and TCRVγnon9 γδ cells infiltrated the pretreated BCC and MEL (on average 2.5 cells/pre-Tt sample), the TCRVγnon9 subset increased more significantly after therapy (on average: 5.3 TCRVγnon9 and 3.6 TCRVγ9 cells/ post-Tt sample, Student p = .04).

Finally, the patient's clinical responses did not correlate to either of i) γδ T and/or T CD8 TIL counts alone or together, ii) TIL counts per milestone, iii) pre-treatment infiltration pattern, or iv) infiltrate reactivity to treatment, taken alone. Thus, beyond confirming that ICB therapy increases the tumorinfiltrating T CD8 and Tex regardless of clinical outcome, [START_REF] Yost | Clonal replacement of tumor-specific T cells following PD-1 blockade[END_REF] our results demonstrated that such treatments affected the γδ T as well as T CD8 TILs.

Discussion

Understanding whether the functional differentiation of tumor-infiltrated γδ T lymphocytes is similar or unrelated to cancer type). d, bubble plot of Tex and Ttrm cell counts per tumor sample, showing cancer type (color key as c), and infiltration pattern. e, representative example of trajectory maps of TILs from the inflamed tumor HNSCC#26, overlaid on the PBMC-derived cells (gray). Color key: TCRVγ9 cells (light blue), Tex TCRVγ9 cells (dark blue), TCRVγnon9 cells (orange), Tex TCRVγnon9 cells (dark red), T CD8 (yellow), Tex T CD8 (black). f, rates of Tex cells per differentiation stage and T cell lineage. Cell counts are indicated. g, violin plots and single cell scores of glycolytic gene enrichment in exhausted and non-exhausted cells from the specified group (Mann-Whitney P values). h, association of exhaustion with tissue residency in the Tcm and Tem stages of T γδ and T CD8 TILs (χ2 P values), absolute counts are specified. Abbreviations: Tn: naïve T, Tcm: central memory T, Tem: effector memory T, Temra: effector memory CD45RA T, Tex: exhausted T, Ttrm: tissue-resident memory. that of the T CD8 compartment is critical to implement γδ T cell-targeted cancer therapies. Here, we provided the most comprehensive cancer-wide atlas of these lymphocytes visualized in their continuous dynamics of differentiation. Thanks to the high resolution of pseudotimed trajectories, our study refined more and subtler transitions than the former four Tnto-Temra stages of differentiation [START_REF] Angelini | FcgammaRIII discriminates between two subsets of Vgamma9Vdelta2 effector cells with different responses and activation pathways[END_REF][START_REF] Caccamo | Differential requirements for antigen or homeostatic cytokines for proliferation and differentiation of human Vgamma9Vdelta2 naive, memory and effector T cell subsets[END_REF] in γδ T lymphocyte and T CD8 lymphocytes. By quantifying each cell's progression along a continuous series of differentiation states, pseudotime represents a transform of true chronological time albeit through non-constant scales. [START_REF] Alpert | Alignment of single-cell trajectories to compare cellular expression dynamics[END_REF][START_REF] Trapnell | The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells[END_REF][START_REF] Van Den Berge | Trajectory-based differential expression analysis for single-cell sequencing data[END_REF] How explicit time is quantitatively converted into T cell pseudotime remains to be determined, but the related sizes of the γδ T and T CD8 trajectories suggests both lineages display the same time transform. Nonetheless, the 'public' trajectory paved here by γδ T lymphocytes from ~150 healthy and diseased individuals (Tables S3,S5) constitutes a reference roadmap of γδ T cell differentiation. Upon direct integration [START_REF] Stuart | Comprehensive integration of single-cell data[END_REF] of newer γδ T cells from any other scRNAseq dataset, each single cell coordinates on this trajectory will define its differentiation hallmarks, allowing automated characterizations based on this resource.

Although well characterized in the mouse, [START_REF] Fiala | From thymus to periphery: molecular basis of effector gammadelta-T cell differentiation[END_REF] the differentiation of peripheral human γδ T cell has remained less investigated over the last decade. A linear sequence of four successive and non-overlapping Tn to Temra stages initially established for TCRVγ9 cells [START_REF] Dieli | Differentiation of effector/memory Vdelta2 T cells and migratory routes in lymph nodes or inflammatory sites[END_REF] was assumed to also apply for TCRVγnon9 cells. This model was validated recently at the single cell level, [START_REF] Pizzolato | Single-cell RNA sequencing unveils the shared and the distinct cytotoxic hallmarks of human TCRVdelta1 and TCRVdelta2 gammadelta T lymphocytes[END_REF] but neither its pseudotime scale nor milestones were known. Here, these features were defined for both subsets, and enlighted several subtle molecular transitions among Tcm and Tem toward stem-like, migration-prone, and exhausted cells, which so far were only characterized in the T CD8 lineage. Hence, most if not all T CD8 differentiation stages also exist in the TCRVγ9 and TCRVγnon9 subsets of γδ T lymphocytes. These subsets indeed have different tropisms at the steady state, with the predominance of TCRVγ9 cells in the blood and TCRVγnon9 cells in extra-lymphoid tissues. [START_REF] Marco | Epithelia use butyrophilin-like molecules to shape organ-specific gammadelta T cell compartments[END_REF] This repartition disappears in tumors, however, where the TCRVγnon9 γδ TILs predominate in all (non-)hematopoietic cancer types. [START_REF] Wu | An innate-like Vδ1+γδ T cell compartment in the human breast is associated with remission in triple-negative breast cancer[END_REF][START_REF] Mikulak | NKp46expressing human gut-resident intraepithelial Vdelta1 T cell subpopulation exhibits high antitumor activity against colorectal cancer[END_REF][START_REF] Foord | Characterization of ascites-and tumor-infiltrating gammadelta T cells reveals distinct repertoires and a beneficial role in ovarian cancer[END_REF] Viral infections are major drivers of peripheral differentiation of T CD8 and γδ T lymphocytes. Thus, the locally phased differentiation of both γδ T and T CD8 lymphocytes strongly suggests a functional convergence toward effective antitumor immunity. [START_REF] Tosolini | Assessment of tumorinfiltrating TCRVγ9Vδ2γδlymphocyte abundance by deconvolution of human cancers microarrays[END_REF][START_REF] Wu | An innate-like Vδ1+γδ T cell compartment in the human breast is associated with remission in triple-negative breast cancer[END_REF][START_REF] Savas | Single-cell profiling of breast cancer T cells reveals a tissue-resident memory subset associated with improved prognosis[END_REF] Substantiating this view, the highest correlation of γδ T and T CD8 TILs cell counts and differentiation observed in virus-positive patients is consistent with the repertoire of T CD8 TILs specificities prominently directed to viral antigens. [START_REF] Simoni | Bystander CD8(+) T cells are abundant and phenotypically distinct in human tumour infiltrates[END_REF] It is tempting to speculate that this functional convergence scenario underlies the better outcome of HPV + over HPV -HNSCC patients, [START_REF] Hewavisenti | CD103+ tumor-resident CD8+ T cell numbers underlie improved patient survival in oropharyngeal squamous cell carcinoma[END_REF] the Sars-CoV2induced remission in EBV + HL, [START_REF] Challenor | SARS-CoV-2-induced remission of Hodgkin lymphoma[END_REF] and the generally lower risk of cancer in CMV + than CMV -kidney transplant recipients. [START_REF] Couzi | Cytomegalovirus-induced gammadelta T cells associate with reduced cancer risk after kidney transplantation[END_REF] Since previous transcriptomic deconvolutions of human cancers showed a dissociation of γδ T and αβ T TILs, [START_REF] Tosolini | Assessment of tumorinfiltrating TCRVγ9Vδ2γδlymphocyte abundance by deconvolution of human cancers microarrays[END_REF][START_REF] Meraviglia | Distinctive features of tumor-infiltrating gammadelta T lymphocytes in human colorectal cancer[END_REF][START_REF] Gentles | The prognostic landscape of genes and infiltrating immune cells across human cancers[END_REF] the tight association of γδ T and T CD8 TILS unveiled here was unexpected. This discrepancy reflects not only the better resolution of scRNAseq transcriptomics than artificial intelligencedriven inferences from bulk microarray deconvolutions, but also the difference between T CD8 alone and the total αβ T compartment. Nevertheless, the coordination of γδ T and T CD8 TILs reported here might also reflect a broader orchestration, possibly involving a BTN3A1-mediated mechanism. [START_REF] Payne | BTN3A1 governs antitumor responses by coordinating alphabeta and gammadelta T cells[END_REF] Furthermore, since T cell differentiation ensues antigenic activation, the advanced differentiation of TCRVγ9 and TCRVγnon9 γδ TILs questions the nature of their driving stimuli in tumor. Whether these are metabolites like bacterial phosphoantigens, [START_REF] Fournie | Stimulation of gamma delta T cells by phosphoantigens[END_REF] viral ligands, [START_REF] Dechanet | Implication of gammadelta T cells in the human immune response to cytomegalovirus[END_REF] stress-induced antigens, [START_REF] Marlin | Sensing of cell stress by human gammadelta TCR-dependent recognition of annexin A2[END_REF] or even HLA-restricted mutational peptide neoantigens remains to be determined more systematically.

The correlated exhaustions of intratumoral γδ T and T CD8 lymphocytes reflect a sustained activation and microenvironmental pressure in both cases. As a consequence of this coherence, the tumor infiltration by both cell types was affected by checkpoint blockade therapies. Although these are critical determinants of outcome, neither T CD8 TILs pre-or post-treatment [START_REF] Tumeh | PD-1 blockade induces responses by inhibiting adaptive immune resistance[END_REF] nor their exhaustion status [START_REF] Daud | Tumor immune profiling predicts response to anti-PD-1 therapy in human melanoma[END_REF][START_REF] Riaz | Tumor and microenvironment evolution during immunotherapy with nivolumab[END_REF][START_REF] Prat | Immunerelated gene expression profiling after PD-1 blockade in non-small cell lung carcinoma, head and neck squamous cell carcinoma, and melanoma[END_REF] suffice to robustly predict the clinical responses of cancer patients. Hence, it is not surprising that the treatment-induced changes of γδ T and T CD8 TIL cell counts, alone or together, did not correlate either. In addition to immunotherapy-induced changes in the functional state of TILs, [START_REF] Efremova | Targeting immune checkpoints potentiates immunoediting and changes the dynamics of tumor evolution[END_REF] the anatomical location of metastases, and tumor mutation of the IFNγ signaling or antigen-presentation pathways [START_REF] Patel | Identification of essential genes for cancer immunotherapy[END_REF] are among the other determinants of therapeutic responses. Nevertheless, our recurrent detection of post-treatment TCRVγ9 γδ TILs in tumors without such cells before treatment suggests that ICB therapy harnesses a circulating γδ T cell contingent to the cancer immunity cycle. [START_REF] Chen | Elements of cancer immunity and the cancer-immune set point[END_REF] Since therapeutic antibodies refuel tumor lesions with fresh, functional, and expanded lymphocytes from peripheral blood, [START_REF] Wu | Peripheral T cell expansion predicts tumour infiltration and clinical response[END_REF] it is conceivable that the recruitment of any circulating subset of expanded γδ T lymphocytes allows intratumoral collaborations of synchronous γδ T and T CD8 effectors. [START_REF] Savas | Single-cell profiling of breast cancer T cells reveals a tissue-resident memory subset associated with improved prognosis[END_REF] Overall, our study builds a fundamental framework to characterize γδ T lymphocyte differentiation in human cancers, and represents an invaluable resource for the emerging cancer immunotherapies targeting γδ T lymphocytes.

versus CD39 (right) in γδ T and T CD8 TILs from the specified patients. Right: composite of results for rates of PD-1 hi TIM-3 + and PD-1 hi CD39 + cells among the γδ T and in T CD8 TILs (shown are results from three independent experiments from the five biopsies). c, Pearson correlation between the specified cells (same patients as above). d, e, single cell suspensions from TC1-engrafted murine HNSCC tumors (n = 12) were stained ex vivo for murine CD3, CD4, CD8, pan-TCRγδ, CD45.2, PD-1, TIGIT, CTLA-4 and TIM-3 and analyzed by flow cytometry. d, PD-1 versus TIM-3 expression in γδ T (top) and T CD8 (bottom) TILs from one mouse (representative result). e, proportions of subpopulations expressing the 16 possible combinations of PD-1, TIGIT, CTLA-4 and TIM-3 immune checkpoints in γδ T (top) and T CD8 (bottom) TIL. The five higher median proportions of TIL subpopulations are represented for γδ T (blue) and T CD8 TIL (orange). Composite of results means and sd. Abbreviations: Tn: naïve T, Tcm: central memory T, Tem: effector memory T, Temra: effector memory CD45RA T, Tex: exhausted T, Ttrm: tissue-resident memory.

Materials and Methods

Clinical patient sample collection

The SMZL samples were collected at the Institut Universitaire du Cancer de Toulouse-Oncopole (IUCT-O) and processed at the CRB Cancer from Hôpitaux de Toulouse, following the standard ethical procedures (Declaration of Helsinki), after obtaining written informed consent from all patients and upon approval by the Institutional Review Board (nDC-2009-989). The carcinoma samples were collected from patients with head and neck squamous cell carcinoma, ovarian cancer, and cervical cancer at the time of surgery at the Institut Universitaire du Cancer de Toulouse-Oncopole (IUCT-O) in accordance with the Declaration of Helsinki, upon approval by the Institutional Review Board (nDC-2016-2656 and DECIdE protocol, NCT03958240) and signed informed consent.

Mouse models

Female C57BL/6 mice, 7 weeks old, were purchased from Janvier Laboratory (Le Genest-Saint-Isle, France) and then maintained in pathogen-free facilities for experiments and for the specified analyses. Experimental protocols were approved by the Regional Ethic Committee of Toulouse Biological Research Federation (ref C2EA-01, FRBT) and by the French Ministry for Higher Education and Research. The European directive 2010/63/EU was followed for guidelines on animal welfare.

Cells and flow cytometry

Healthy donors PBMC were obtained after Ficoll-Hypaque (GE Healthcare) density centrifugation and cultured in RPMI 1640 medium (Dutscher) supplemented with 10% fetal bovine serum (ThermoFisher Scientific). Human tumor biopsy specimens were weighed, washed with PBS and mechanically disrupted. These preparations were centrifuged on a Ficoll gradient and frozen. When specified in text, T cells were isolated by negative selection using the Human T Lymphocytes Enrichment Set-DM (BD Biosciences) according to the manufacturer's instructions. TC1 cells expressing the HPV 16 E6 and E7 proteins were developed in the laboratory of T.C. Wu (Johns Hopkins University). Resting PBMC control was processed for immunostaining with the specified antibodies or ADT (Biolegend) for 20 min at room temperature. For intracellular stainings, cells were fixed and permeabilized using the transcription factor staining buffer set (eBioscience) before staining with the proper antibodies or reagent. Used antibodies are listed on the Key Resource Table . Flow cytometry acquisitions were performed on a BD LSR II cytometer (BD Biosciences) and data were analyzed with Cytobank (http:// www.cytobank.org).

Single cell RNA sequencing

The Chromium Controller Single-Cell Instrument and Chromium Single Cell 3′ Library & Gel Bead Kit (V3) were used to prepare individually barcoded single-cell RNA-Seq libraries, following manufacturer's instruction (10x Genomics). Briefly, cells were thawed and viability was improved by using Dead Cell Removal Kit (Miltenyi Biotec). Then, cells were washed twice with calcium and magnesiumfree phosphate buffered saline containing 0.04% bovine serum albumin, and re-suspended at a concentration comprised between 700 and 1200 cells/µL (viability>90%) in the same buffer. In the experiment, 10 000 cells of each sample were loaded to obtain a yield comprised between 2000 and 7000 sequenced cells. Reverse transcription was performed on the C1000 Touch Thermal Cycler, and cDNA was amplified for 12 cycles. cDNA quality control and dosage were performed on Fragment Analyzer, and then by PCR. Double-sided size selection of PCR products was performed with SPRIselect Reagent Kit (Beckman Coulter). The sequencing-ready library was cleaned up with SPRIselect beads (Beckman Coulter) and the sequencing was performed on a HiSeq 3000 instrument (Illumina) with a sequencing depth of at least 50,000 reads per cell.

scRNA-seq data pre-processing and QC

The sequencing output bcl2 files were converted to FASTQ format by using cellranger-mkfastq™ algorithm (10x Genomics), and cellranger-count was applied to align to the GRCh38 reference transcriptome and build the (cell, UMI) expression matrix for each sample. Preprocessing, normalization of UMI, and QC were assessed using the R package Seurat. [START_REF] Stuart | Comprehensive integration of single-cell data[END_REF][START_REF] Butler | Integrating single-cell transcriptomic data across different conditions, technologies, and species[END_REF] First, genes that correspond to human leukocyte antigens (HLA) and genes located in chromosome X and Y were removed, and the quality of cells was then determined by the number of detected genes per cell, and the proportion of mitochondrial genes. Additional datasets downloaded from 10x Genomics website, Gene Expression Omnibus (GEO), European Genome-phenome Archive (EGA), and ArrayExpress followed the same preprocessing and QC analyses.

Gene annotations

All gene symbols and protein annotations were from public databases (Uniprot, Genecards, and The Human Protein Atlas for humans) and transcription factors genes were defined as such from the GO annotations.

Integration of PBMC scRNAseq datasets

Here, seven different scRNAseq datasets were assembled to digitally extract control γδ T lymphocytes from PBMC of several healthy individuals. This comprised (n = 8,202) pure γδ T lymphocytes cell-sorted by cell surface TCRVδ1 or TCRVδ2 from the PBMC of n = 3 healthy donors [START_REF] Pizzolato | Single-cell RNA sequencing unveils the shared and the distinct cytotoxic hallmarks of human TCRVdelta1 and TCRVdelta2 gammadelta T lymphocytes[END_REF] and (n = 25,736) PBMC from four other healthy donors (4k +8k from 10XGenomics-3ʹplatform V2, and 10k +1k from 10XGenomics-3ʹplatform V3 downloaded from 10xGenomics resource). After pre-processing and QC discarding of cell doublets and dying cells, the seven datasets were integrated using Seurat 3.0. [START_REF] Stuart | Comprehensive integration of single-cell data[END_REF] Briefly, the datasets were normalized together and variable features from each dataset were identified by FinVariable-Data before identification of anchors correspondences across all the datasets, to yield an integrated PBMC dataset containing 26,313 features across n = 25,252 control cells (n = 17,803 PBMC plus n = 7,449 purified γδ T lymphocytes). Principal component analysis (PCA) was performed on this integrated dataset and used to calculate uniform manifold approximation and projection (UMAP), as well as t-distributed Stochastic Neighbor Embedding (t-SNE) [START_REF] Van Der Maaten | Visualizing high-dimensional data using t-SNE[END_REF] coordinates. When visualized in UMAP space, [START_REF] Mcinnes | Umap: Uniform manifold approximation and projection for dimension reduction[END_REF] its main cell types were annotated by expression of the CD14, CD19, CD3D, and NCAM marker genes (Supplemental Data Figure 1a).

Single cell signatures and scores

Integrated and normalized scRNAseq datasets were then scored for the multigene signatures specified in text (Table S1) which were downloaded from the MSigDB (https:// www.gsea-msigdb.org/gsea/msigdb/collections.jsp) unless specified otherwise. Each single cell score for each geneset was computed as described. [START_REF] Pont | Single-cell signature explorer for comprehensive visualization of single cell signatures across scRNA-seq datasets[END_REF] Briefly, Single-Cell Signature Scorer computes the score of cell C j for gene set GS x as the sum of all UMI for all the GS x genes expressed by C j divided by the sum of all UMI expressed by Cj: Score of cell C j for geneset GS x = (∑(UMI) GSx )/(∑(UMI) Cj ) [START_REF] Dieli | Differentiation of effector/memory Vdelta2 T cells and migratory routes in lymph nodes or inflammatory sites[END_REF] For further visualization, all cell signatures scores are merged with each cell's UMAP coordinates using Single-Cell Signature Merger. [START_REF] Pont | Single-cell signature explorer for comprehensive visualization of single cell signatures across scRNA-seq datasets[END_REF]

Score and Gate of γδ T and T CD8 cells from PBMC

Based on using single-cell signature scores to identify cells, this procedure allows the selective and exhaustive digital isolation from any scRNAseq dataset of any subset of cells specified by characteristic single-cell scores. The γδ T cell detection, which previously relied on a singlegene signature, [START_REF] Pizzolato | Single-cell RNA sequencing unveils the shared and the distinct cytotoxic hallmarks of human TCRVdelta1 and TCRVdelta2 gammadelta T lymphocytes[END_REF] was improved here by implementing a sequence of positive and negative cell selections from several gene signatures, called 'Score and Gate.' The above (n = 25,252) PBMC-extracted cells were scored for the specified gene signatures (Table S1) using Single-Cell Signature Explorer [START_REF] Pont | Single-cell signature explorer for comprehensive visualization of single cell signatures across scRNA-seq datasets[END_REF] and scores were sequentially visualized on scatterplots by Single-Cell Virtual Cytometer. [START_REF] Pont | Single-cell virtual cytometer allows user-friendly and versatile analysis and visualization of multimodal single cell RNAseq datasets[END_REF] The Score and Gate procedure comprises five sequential gating steps, here applied starting from the (n = 17,803) PBMC sub-dataset alone:

1. Positive selection of double negative cells [= T and NK lineages (n = 11,768)] from the scatterplot of 'B cell' and 'Myeloid cell' signatures.

2. Scatterplot of 'B cell' signature against the addition of 'CD8AB' and 'TCRγδ' signatures to positively select the double positive cells [= non-CD4 T cells (n = 6,162)].

3. Scatterplot of 'CD8ab' signature against 'TCRα constant gene TRAC' signature to digitally extract doublepositive cells [= αβ T CD8 lymphocytes (n = 3,680)], on the one hand, and negative cells [= γδ T (n = 2,161) plus CD4CD8-double negative αβ T cells], on the other hand.

4. Scatterplot of the addition of 'CD3 complex' and 'TCRα constant gene TRAC' signatures against the addition of 'CD3 complex' and 'TCRγδ' signatures to extract the 'CD3 complex'-'TCRγδ' double-positive cells [=γδ T lymphocytes (n = 580, 3% of initial PBMC)].

TCRVγ9 and TCRVγnon9 cell subtyping by compensated scores

Since γδ T lymphocytes express either the TRGC1-encoded TCRVγ9 or the TRGC2-encoded TCRVγnon9 in a mutually exclusive fashion, the extracted γδ T lymphocytes were classified on 'TRDC, TRGC1ʹ positive cells [=TCRVγ9 cells (n = 312)] or 'TRDC, TRGC2ʹ positive cells [=TCRVγ non 9 cells (n = 268)] using a compensated version of these two signatures. Compensated scores are obtained by multiplying the score of the geneset of interest, here GSS TCRVγ9 , by its difference to its complementary geneset (GSS TCRVγ9 minus GSS TCRVγnon9 ): [START_REF] Dieli | Differentiation of effector/memory Vdelta2 T cells and migratory routes in lymph nodes or inflammatory sites[END_REF] Scatterplots of compensated scores for the two complementary genesets GSS TCRVγ9 and GSS TCRVγnon9 unambiguously identify the TCRVγ9 or TCRVγnon9 subset of each γδ T lymphocyte (Supplemental Data Figure 1b).

GSS TCRVγ9 comp = GSS TCRVγ9 * (GSS TCRVγ9 -GSS TCRVγnon9 )

Validation of the Score and Gate method of digital extraction

The method was validated on two different scRNAseq datasets: our formely depicted dataset of purified γδ T lymphocytes from three healthy donors [START_REF] Pizzolato | Single-cell RNA sequencing unveils the shared and the distinct cytotoxic hallmarks of human TCRVdelta1 and TCRVdelta2 gammadelta T lymphocytes[END_REF] and a CITEseq dataset of PBMC from another healthy donor. [START_REF] Pont | Single-cell virtual cytometer allows user-friendly and versatile analysis and visualization of multimodal single cell RNAseq datasets[END_REF] When applied to the (n = 7,449) purified γδ T cell sub-dataset, this same score and gate strategy extracted (n = 5,542) γδ T cells (74% sensitivity). Given their already known TCRVδ1/ TCRVδ2 subtype, the TCRVγ9/non9 subdivision delineated (n = 2,328) TCRVγ9Vδ2, (n = 1,993) TCRVγnon9Vδ1, (n = 817) TCRVγnon9Vδ2, and (n = 404) TCRVγ9Vδ1 cells. Notwithstanding the equal rates of annotated TCRVδ1 and TCRVδ2 cells deliberately composing this sub-dataset, [START_REF] Pizzolato | Single-cell RNA sequencing unveils the shared and the distinct cytotoxic hallmarks of human TCRVdelta1 and TCRVdelta2 gammadelta T lymphocytes[END_REF] such results confirmed the prominence of TCRVγ9Vδ2 cells over all the other TCR subsets of blood γδ T, and the preferential pairings of TCRVγ9 (85%) with TCRVδ2, and of TCRVγnon9 (71%) with TCRVγ1 (Supplemental Data Figure 1c). In the UMAP, (n = 5,537 out of the 5,542) extracted γδ T cells are embedded in the T cell clusters (Supplemental Data Figure 1d) showing the γδ T cell detection specificity by score and gate had improved (99.9 vs 95%) over earlier method, [START_REF] Pizzolato | Single-cell RNA sequencing unveils the shared and the distinct cytotoxic hallmarks of human TCRVdelta1 and TCRVdelta2 gammadelta T lymphocytes[END_REF] and successfully extracted those γδ T cells embedded in CD4 and CD8 T cells clusters. Since CITEseq datasets provide both transcriptomes and immune phenotypes of single cells, [START_REF] Stoeckius | Simultaneous epitope and transcriptome measurement in single cells[END_REF] the same score and gate sequence was also applied to our CITEseq dataset of 3k PBMC from a healthy donor (GSE144434). [START_REF] Pont | Single-cell virtual cytometer allows user-friendly and versatile analysis and visualization of multimodal single cell RNAseq datasets[END_REF] This extracted (n = 137) γδ T cells (4.7% of PBMC, encompassing n = 81 TCRVγ9 plus n = 56 TCRVγnon9 cells) of cell surface CD3 + CD4 -CD8 - phenotype, and embedded with T CD8 and NK cells in UMAP (Supplemental Data Figure 7), as expected. These results validated the score and gate method, and the lymphocytes extracted from the CITEseq dataset were integrated into the control datasets. The integrated (cell x gene) matrices of PBMC-derived γδ T and T CD8 single cells are provided in Tables S3, S4, respectively, together with their annotations for pseudotime and differentiation hallmarks in Tables S5, S6.

Extraction of TILs from scRNAseq datasets of tumor samples

The scRNA-seq of tumor lesions from three patients with splenic marginal zone lymphoma (SMZL) were produced for this study, and additional datasets of tumors from other cancers were downloaded from public resources (GEO, EGA, and AE). These encompassed tumor samples from cutaneous T cell lymphoma (CTCL), [START_REF] Gaydosik | Single-cell lymphocyte heterogeneity in advanced cutaneous T-cell lymphoma skin tumors[END_REF] Hodgkin lymphoma (HL), [START_REF] Aoki | Single-cell transcriptome analysis reveals disease-defining T-cell subsets in the tumor microenvironment of classic Hodgkin lymphoma[END_REF] lung adenocarcinoma (LUAD) [START_REF] Lambrechts | Phenotype molding of stromal cells in the lung tumor microenvironment[END_REF][START_REF] Kim | Single-cell RNA sequencing demonstrates the molecular and cellular reprogramming of metastatic lung adenocarcinoma[END_REF] , breast carcinoma (BRCA), [START_REF] Azizi | Singlecell map of diverse immune phenotypes in the breast tumor microenvironment[END_REF] triple negative breast cancers (TNBC), [START_REF] Savas | Single-cell profiling of breast cancer T cells reveals a tissue-resident memory subset associated with improved prognosis[END_REF] head and neck squamous cell carcinoma (HNSCC), [START_REF] Cillo | Immune landscape of viral-and carcinogen-driven head and neck cancer[END_REF] hepatocellular carcinoma (HCC), [START_REF] Zheng | Landscape of infiltrating T cells in liver cancer revealed by single-cell sequencing[END_REF] basal cell carcinoma (BCC), [START_REF] Yost | Clonal replacement of tumor-specific T cells following PD-1 blockade[END_REF] and melanoma (MEL) [START_REF] Sade-Feldman | Defining T cell states associated with response to checkpoint immunotherapy in melanoma[END_REF] (Table S2). The same score and gate procedure was applied to the pre-processed and quality controlled datasets for digital extraction of (n = 4,648) γδ T and (n = 45,174) T CD8 cells, referred to as TILs. From the CTCL samples, the malignant T lymphocytes were filtered out from TILs by gating out those cells overexpressing the marker gene TOX and the CTCL-specific 17 gene signature. [START_REF] Gaydosik | Single-cell lymphocyte heterogeneity in advanced cutaneous T-cell lymphoma skin tumors[END_REF] The integrated (cell and gene) matrices of both PBMC-and TILs-derived γδ T and T CD8 single cells are, respectively, provided in Tables S3, S4

with their respective annotations in Tables S5,S6.

Trajectory Inference After the digital extraction of γδ T and T CD8 lymphocytes from normal and cancer samples, both control and TILs were integrated, independently for each lineage (γδ T and T CD8), by using Seurat (R package), [START_REF] Stuart | Comprehensive integration of single-cell data[END_REF] as described above. Next, the γδ T and T CD8 integrated matrices were independently used for maturation trajectory inference using dynverse (collection of R packages), which provides more than 45 different methods for trajectory inference (TI). The TI was calculated following the recommended workflow. [START_REF] Saelens | A comparison of single-cell trajectory inference methods[END_REF] Briefly, the normalized (log2) of raw expression data, and UMAP coordinates were used as an input. The most suitable and rapid method for our datasets determined by applying 'guidelines' function (dynverse on R) was Minimum Spanning Tree (MST) to infer the maturation trajectories of the γδ T and T CD8, independently ('ti_mst' function, dynverse on R), with default settings. Additional signatures ('Tn' and 'Teff') were scored, as described above, and together with the obtained MST attributes: dimension reduction coordinates 1 and 2 (MST1 and MST2), milestones and pseudotime, were used to define the root of the trajectory ('add_root' function) of the γδ T and T CD8 lymphocytes. The MST attributes allow multiple visualizations, so the data were visualized through pseudotimed trajectories.

Pseudotimed trajectories

The pseudotimed trajectory is a bidimensional representation of the trajectory inferred for any single cell (p) in which (x,y,z) are the MST coordinates for pseudotime, MST1, and MST2, respectively. The dim_proj method takes each cell's (x,y,z) from this 3-dimension space (X: pseudotime, Y: MST1, Z: MST2) and projects them in a 2-dimensional plane where pseudotime is one of these dimensions (XoY'). Between the three orthonormal dimensions (XYZ) and the projection plane (XoY'), there is a rotation angle (ZoY') called β (Supplemental Data Figure 8a). The new projected coordinates (x,y') of the cell p are calculated as follows:

o'p 2 = y 2 + z 2 o'y' = ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi

ðy þ y min Þ 2 þ ðz þ z min Þ 2 q * cos(arctan ((y þ y min Þ= z þ z min ð Þ)+α) (5)
The resulting trajectory using only cell's coordinates (x, y') is called pseudotimed trajectory. Of note, although the trajectory itself is independent of the rotation angle β, selecting this latter allows us to optimize visualization of the same trajectory and its eventual branches (Supplemental Data Figure 8b).

Heatmaps

Heatmaps of gene expression level (UMI mean) per milestone were generated with color scale based on the z-score distribution per row, and shown from row min to row max. Heatmaps of multigene signature score per milestone (milestone mean) were generated by computing scores with Single-Cell Signature Explorer [START_REF] Pont | Single-cell signature explorer for comprehensive visualization of single cell signatures across scRNA-seq datasets[END_REF] and shown by color scale from the z-score distribution per row. Heatmaps of cell counts per milestone were generated as fractions of the row total.

Classifications

The cell classifications Ttrm and Tex were performed by 'at least one binary' as follows. The specified single cells were scored for several reference gene signatures (Table S1). For Tex classification, the reference Tex signatures were five published and partially overlapping Tex genesets. [START_REF] Khan | TOX transcriptionally and epigenetically programs CD8(+) T cell exhaustion[END_REF][START_REF] Balanca | Dual relief of T-lymphocyte proliferation and effector function underlies response to PD-1 blockade in epithelial malignancies[END_REF][START_REF] Chihara | Induction and transcriptional regulation of the co-inhibitory gene module in T cells[END_REF][START_REF] Alfei | TOX reinforces the phenotype and longevity of exhausted T cells in chronic viral infection[END_REF][START_REF] Tosolini | Large-scale microarray profiling reveals four stages of immune escape in non-Hodgkin lymphomas[END_REF] For each signature, a cutoff was defined as the maximal score of the (n = 3,680) control γδ T CD8 lymphocytes extracted from healthy adults PBMC. This threshold defined the cell's binary (1: cell score > threshold, 0: otherwise). The Tex cutoffs were: 3.9 for 'Chihara_IL27_Coinhib_module'; 61 0.52 for 'Alfei_d20_tox'; 62 0.22 for 'Khan_Tox_OverExpressed_genes'; 16 0.5 for 'Tosolini_ NHL_IEGS33ʹ; [START_REF] Tosolini | Large-scale microarray profiling reveals four stages of immune escape in non-Hodgkin lymphomas[END_REF] and 0.16 for 'Balanca_QP_genes'. [START_REF] Balanca | Dual relief of T-lymphocyte proliferation and effector function underlies response to PD-1 blockade in epithelial malignancies[END_REF] For each TIL, the five Tex binaries were summed, and any TIL with a non-zero sum of binaries was classified Tex. The Ttrm classification was applied likewise using six reference and partially overlapping Ttrm signatures from the literature. [START_REF] Wu | Peripheral T cell expansion predicts tumour infiltration and clinical response[END_REF][START_REF] Kumar | Human tissue-resident memory T cells are defined by core transcriptional and functional signatures in lymphoid and mucosal sites[END_REF] The binarizing cutoffs established with the control T CD8 lymphocytes were: 0.46 for 'Kumar_13g_Ttrm'; 0.3 for 'Kumar_3g_Ttrm'; [START_REF] Kumar | Human tissue-resident memory T cells are defined by core transcriptional and functional signatures in lymphoid and mucosal sites[END_REF] 3.72 for 'Wu_Tcellcluster4.1_trm'; 2.6 for 'Wu_Tcellcluster8.3_trm'; 2.8 for 'Wu_Tcellcluster8.3b_trm'; 4.2 for 'Wu_Tcellcluster8.3c_trm'. [START_REF] Wu | Peripheral T cell expansion predicts tumour infiltration and clinical response[END_REF] These Ttrm binaries were summed, and any TIL with a non-zero sum of Ttrm binaries was classified Ttrm. The cytolytic T cell infiltrates of tumors were arbitrarily classified by total cytolytic (T CD8 + γδ T) TIL counts as deserts (D, <100 TILs), classic (C, 100-1000 TILs), or inflamed (I, >1000 TILs). The Tex classification from the flow cytometry data was based on both median (PD-1 hi TIM-3 + )>8% and median (PD-1 hi CD39 + )>12%. The reactivity of tumor infiltrates (TI) from ICBtreated MEL and BCC patients was classified as TI-R (reactive)/TI-NR (non-reactive) based on the difference of total TIL counts (post-pre-ICB therapy) by TI-R: ∆>0, else TI-NR.

Figure 1 .

 1 Figure 1. Circulating γδ T are more differentiated than T CD8 lymphocytes in most healthy individual's blood. a, pseudotimed trajectory of PBMC-derived γδ T cells. Cell color key: pseudotime-ordered milestones. b, mean of gene expression per milestone. Genes shown blue and red: naïve (Tn) and effector (Teff) signatures, respectively. c, Single cell scores for Tn and Teff signatures across the trajectory. d, dendrogram of milestones (centroids, same color key as a) overlaid to the pseudotimed trajectory showing all γδ T cells (gray dots), marker genes, trajectory (dark line), and differentiation direction (lines). e-h, same as a-d for T CD8 cells extracted from the same PBMC. i, rates of γδ T and T CD8 cells at each differentiation stage. Abbreviations: Tn: naïve T, Tcm: central memory T, Teff: T effector, Tcm-mig: migration-prone central memory T, Tem: effector memory T, Tem-act: activated effector memory T, Tem-cyto: cytotoxic effector memory T, Temra: effector memory CD45RA T, Tscm: stem-like central memory T, Tpex: progenitor exhausted-like T.

Figure 2 .

 2 Figure 2. Memory γδ T and CD8 T cells are enriched in 136 human tumors from 11 cancer types. a, pseudotimed trajectories of TILs (red) of the (n = 4,648) γδ T and (n = 45,174) T CD8 cells featuring their dendrogram of milestones, overlaid to their PBMC-derived counterpart (green) shown on their respective trajectory. b, Rates of γδ T (blue) and T CD8 (orange) TILs with the specified differentiation stage in representative tumor samples. c, TIL cell counts and rates in each tumor sample (colored by

Figure 3 .

 3 Figure 3. The ex vivo phenotype of tumor-infiltrating γδ T and CD8 T lymphocytes consists of similarly differentiated cells a, right: Cell surface CCR7 and CD45RA phenotype of the γδ T and T CD8 TILs from a freshly biopsied HNSCC tumor, and differentiation stages in this tumor. Left: Rates of Tn, Tcm, Tem, and Temra cells in γδ T and T CD8 lineages from the five freshly biopsied tumors (means from two independent experiments). b, cell surface phenotype of PD-1 versus TIM-3 (left) and PD-1

Figure 4 .

 4 Figure 4. Both γδ T and CD8 T TILs respond to anti-PD1 therapy in patients with cutaneous cancers. a, summary of characteristics and TIL counts extracted from the specified BCC and MEL samples. Inflamed (dark red), classic (pink), desert infiltrates (white). b, representative examples of inflamed (BCC#009), classic (BCC#001), and desert (MEL#008) pre-treatment infiltrate patterns. The T TILs extracted from each specified tumor are overlaid on their respective pseudotimed trajectory (gray: control blood cells of the corresponding lineage; TCRVγ9 cells (light blue), Tex TCRVγ9 cells (dark blue), TCRVγnon9 cells (orange), Tex TCRVγnon9 cells (dark red), T CD8 (yellow), Tex T CD8 (black). Comparison of BCC samples before and after ICB treatment for: c, correlation of T CD8 and γδ T cell counts, d, distribution of differentiation stages in T CD8 and γδ T TILs. Abbreviations: Tn: naïve T, Tcm: central memory T, Teff: T effector, Tcm-mig: migration-prone central memory T, Tem: effector memory T, Tem-act: activated effector memory T, Tem-cyto: cytotoxic effector memory T, Temra: effector memory CD45RA T, Tscm: stem-like central memory T, Tpex: progenitor exhaustedlike T.

þ z 2 p

 2 (1) cos(α+β) = o'y'/o'p(2) So the new Y coordinate y' of p projection in the plane is:o'y' = ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiðy 2 Þ * cos(α+β) (3) o'y' = ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ðy 2 þ z 2 pÞ * cos(arctan(ðy=z)+β) (4) To avoid negative values, the whole trajectory was translated by adding the minimal coordinates (y min , z min ) of the trajectory cells to the initial (y, z) coordinates:

e1939518-8J.-P. CERAPIO ET AL.

Acknowledgments

We acknowledge F. Dieli and S. Coffelt for suggestions and critical comments on the manuscript, C Steidl for facilitating access to HL datasets, and CRCT colleagues for helpful discussions.

Funding

This work was supported in part by institutional grants from the Institut National de la Santé et de la Recherche Médicale (INSERM), the Université Toulouse III-Paul Sabatier, the Centre National de la Recherche Scientifique (CNRS); Laboratoire d'Excellence Toulouse Cancer (TOUCAN-2) under grant contract ANR11-LABX, the Fondation ARC under grant contract PGA1-RF2019-0208691; Institut Universitaire du Cancer de Toulouse-Oncopole under contract CIEL; Institut Carnot Lymphome under cointract CALYM. This work was granted access to the HPC resources of CALMIP supercomputing center under allocation 2020-T19001. We are also grateful to the Genotoul bioinformatics platform (Bioinfo Genotoul, Toulouse Midi-Pyrenees) for providing computing resources;Fondation ARC pour la Recherche sur le Cancer [PGA1 RF20190208691];Institut Carnot CALYM [RAC21009BA]; LabEx Toulouse Cancer [ANR11-LABX].

Data visualization

Single-cell Multilayer Viewer, was developed as a serverless software for the visualization of up to 5 layers of information. It is available on GitHub repository: https://github.com/ MarionPerrier/scMLV. Single-cell Multilayer Viewer allows us to merge and visualize up to 5 layers of quantitative and qualitative variables, alone and in combination. Code and software availability Single-cell Multilayer Viewer (scMLV) developed for this study is available on GitHub repository: https://github. com/MarionPerrier/scMLV Single-Cell Signature Scorer and Single-Cell Virtual Cytometer are available at: https://sites.google.com/site/fredsoftwares/products. Seurat v3.2.2 is available from the Comprehensive R Archive Network (CRAN) and further details in the installation could be found on https://satijalab. org/seurat/. The collection of R packages (dyno v0.1.2 and tidyverse v1.3.0) is available from Github via devtools, with details for the installation from: https://dynverse.org/.

Abbreviations

BCC: basal cell carcinoma, BRCA:Breast carcinoma, CC: cervical cancer, CITEseq: cellular indexing of transcriptomes and epitopes by sequencing, CMV: cytomegalovirus, CTCL: cutaneous T cell lymphoma, EBV: Epstein-Barr virus, HCC: hepatocellular carcinoma, HL: Hodgkin's lymphoma, HNSCC: head and neck squamous cell carcinoma, HPV: human papilloma virus, ICB: immune checkpoint blockade, LUAD: lung adenocarcinoma, MEL: melanoma, OVCAR: ovarian carcinoma, scRNAseq: single-cell RNA sequencing, SMZL: splenic marginal zone lymphoma, Tcm: central memory T, Teff: T effector, Tem: effector memory T, Temra: effector memory CD45RA T, Tpex: progenitor exhausted-like, Tex: exhausted T, TIL: tumor-infiltrated lymphocyte, Tn: naïve T, TNBC: triple-negative breast carcinoma, TnEnCD8: naive-enriched T CD8, Ttrm: tissueresident memory T.

Competing interests statement

Authors declare no relevant conflict of interest.

Data and code availability

Lead contact Further information and requests for resources and reagents should be directed to and will be fulfilled upon reasonable requests by the lead contact, Jean-Jacques Fournie Email jean-jacques.fournie@inserm.fr.

Materials availability This study did not generate new unique reagents.

Data availability The raw data produced in this study are available at Gene Expression Omnibus (GEO) as dataset GSE175785. Published scRNAseq datasets used in this study are available at GEO with the accession numbers: GSE110686, GSE114724, GSE120575, GSE123813, GSE128531, GSE139324, GSE140228; at the European Genomephenome Archive (EGA) with the accession number: EGAS00001004085; and at the Array Express Archive (AE) with the accession number: E-MTAB-6149. The integrated (cell, UMI) matrices of extracted γδ T and T CD8 are provided as Tables S3, S4 and the corresponding cell annotations are available in Tables S5,S6.

Author contributions

Source data

The raw data produced in this study are available at Gene Expression Omnibus (GEO) with accession number GSE175785, and GSE. The integrated (cell, UMI) matrices of all extracted γδ T and T CD8 are provided as Tables S3, S4, and their corresponding cell annotations are in Tables S5, S6.