Expanding the Toolbox for Bicelle-Forming Surfactant–Lipid Mixtures

Bicelles are disk-shaped models of cellular membranes used to study lipid–protein interactions, as well as for structural and functional studies on transmembrane proteins. One challenge for the incorporation of transmembrane proteins in bicelles is the limited range of detergent and lipid combinations available for the successful reconstitution of proteins in model membranes. This is important, as the function and stability of transmembrane proteins are very closely linked to the detergents used for their purification and to the lipids that the proteins are embedded in. Here, we expand the toolkit of lipid and detergent combinations that allow the formation of stable bicelles. We use a combination of dynamic light scattering, small-angle X-ray scattering and cryogenic electron microscopy to perform a systematic sample characterization, thus providing a set of conditions under which bicelles can be successfully formed.

Introduction
Bicelles are nano-sized discoidal-like particles composed of long-chain phospholipids, which form a lipid bilayer core, and detergents or short-chain phospholipids, which form the flanking rim. These structures spontaneously assemble by mixing phospholipids and a detergent at a given lipid:detergent molar ratio (q), which determines the particle size, and their properties vary depending on experimental conditions, such as lipid and detergent type and concentration, buffer ionic strength, and temperature [1][2][3]. Their high compositional tunability, along with stability against several freeze-thawing cycles, makes bicelles popular as models of cellular membranes for the study of both transmembrane protein structure and functioning and their intermolecular interactions in a native-like environment [4,5].
The size of bicelles can be tuned by altering the lipid:detergent molar ratio (q), with the size increasing at higher q values. At q > 1, bicelles become more elongated and such elongated structures make bicelles suitable for nuclear magnetic resonance (NMR) studies due to the alignment of bicelles along the magnetic field [1,6].
Despite the wide range of applications for bicelles as models of biomembranes, detailed structural studies exist only for a handful of compositions [2,17]. The standard bicelle model is composed of a disk with three layers (phospholipid head/tail/head) and a rim (detergent head), but more detailed models have been produced to mimic different regions of the lipid core with varying molecular packings [2]. In this work, we applied DLS, SAXS, and cryo-EM to characterize bicelles developed using either saturated (DMPC/DMPG) or unsaturated 2-oleoyl-1-palmitoyl-sn-glyecro-3-phosphocholine (POPC) and 1-palmitoyl-2-oleoyl-sn-glycero-3-phospho-(1′-rac-glycerol) (POPG) phospholipids combined with three different detergents (DHPC, CHAPS, and sodium cholate). We identify the conditions suitable for controlling the particle size and shape of bicelles composed of the components shown in Figure 1.

Results
The aim of this study was to systematically characterize bicelles composed of different combinations of lipids and detergents to determine the conditions under which stable bicelles can be formed and potentially used for the in-solution study of transmembrane proteins. In this study three different detergents were investigated: CHAPS and sodium

Results
The aim of this study was to systematically characterize bicelles composed of different combinations of lipids and detergents to determine the conditions under which stable bicelles can be formed and potentially used for the in-solution study of transmembrane proteins. In this study three different detergents were investigated: CHAPS and sodium cholate, which were chosen as they have been widely used for the purification of membrane proteins, and DHPC, which has been extensively used to form bicelles in the past.
For each of the three detergents, bicelles were prepared with four lipid compositions: (i) saturated zwitterionic (DMPC), (ii) saturated zwitterionic/anionic (DMPC:DMPG 80:20), (iii) unsaturated zwitterionic (POPC), and (iv) unsaturated zwitterionic/anionic (POPC:POPG 80:20) (Figure 1). Each detergent-lipid combination was prepared at seven different q ratios: 0.5, 1, 1.5, 2, 2.5, 3, and 5 by resuspending lipid/detergent films in a buffer containing the detergent of choice at a concentration slightly above the critical micelle concentration (CMC) and subjecting them to four freeze-thawing cycles to induce bicelle formation prior to structural analysis (see the methods section for details). Figure 2 shows a summary of the hydrodynamic radii (R H ) obtained from the DLS experiments and the radii of gyration (R g ) obtained from the SAXS measurements of the samples. By exploring these phospholipids, we aimed to highlight the effects of the different natures of the tail and headgroup regions on the final assemblies obtained, a first step toward creating more representative models of mammalian and bacterial membranes. cholate, which were chosen as they have been widely used for the purification of membrane proteins, and DHPC, which has been extensively used to form bicelles in the past. For each of the three detergents, bicelles were prepared with four lipid compositions: (i) saturated zwitterionic (DMPC), (ii) saturated zwitterionic/anionic (DMPC:DMPG 80:20), (iii) unsaturated zwitterionic (POPC), and (iv) unsaturated zwitterionic/anionic (POPC:POPG 80:20) (Figure 1). Each detergent-lipid combination was prepared at seven different q ratios: 0.5, 1, 1.5, 2, 2.5, 3, and 5 by resuspending lipid/detergent films in a buffer containing the detergent of choice at a concentration slightly above the critical micelle concentration (CMC) and subjecting them to four freeze-thawing cycles to induce bicelle formation prior to structural analysis (see the methods section for details). Figure  2 shows a summary of the hydrodynamic radii (RH) obtained from the DLS experiments and the radii of gyration (Rg) obtained from the SAXS measurements of the samples. By exploring these phospholipids, we aimed to highlight the effects of the different natures of the tail and headgroup regions on the final assemblies obtained, a first step toward creating more representative models of mammalian and bacterial membranes. DLS measurements showed that most phospholipid/detergents mixtures investigated formed a monodisperse population of particles ( Figure S1) with an increasingly large RH that grew as a function of q. As q was inversely proportional to the molar fraction of detergent in the mixture, higher q values were expected to promote the formation of larger lipid aggregates, in line with previous results [17]. The particle RH followed the trend DHPC > CHAPS > cholate for all phospholipid mixtures at nearly all q values investigated (except for the POPC/POPG-CHAPS mixtures at q 2.5 and 3, which exhibited a larger RH than the DHPC mixtures). For SAXS-based Rg, no estimation could be made for DLS measurements showed that most phospholipid/detergents mixtures investigated formed a monodisperse population of particles ( Figure S1) with an increasingly large R H that grew as a function of q. As q was inversely proportional to the molar fraction of detergent in the mixture, higher q values were expected to promote the formation of larger lipid aggregates, in line with previous results [17]. The particle R H followed the trend DHPC > CHAPS > cholate for all phospholipid mixtures at nearly all q values investigated (except for the POPC/POPG-CHAPS mixtures at q 2.5 and 3, which exhibited a larger R H than the DHPC mixtures). For SAXS-based R g , no estimation could be made for the conditions where large aggregates were present due to the lack of the Guinier regime. The q-dependent growth was considerably less marked for the sodium cholate-containing samples in comparison with DHPC and CHAPS for both saturated and unsaturated lipids. In all sodium cholate mixtures, the R H values increased from a minimum of~15 Å at q = 0.5 to a maximum of~25 Å at q = 5. These values were too small to indicate the presence of a phospholipid bilayer structure; therefore, it was not likely that bicelles formed under these conditions, even at the highest q value tested ( Figure 2, upper panels).
The DLS data did not provide any information regarding the internal structure of the observed particles besides the average R H . Therefore, it is not possible, based on these data alone, to determine what type of particles formed in the solution and whether they were bicelles or mixed micelles [18]. The R H values and the q-dependent behavior of the particle sizes obtained from the DLS data suggest that, with sodium cholate, bicelles do not form in the presence of neither saturated nor unsaturated phospholipids. For CHAPS and DHPC, however, small discoidal bicelles could be representative of the structures formed with saturated lipids, whereas, in the presence of unsaturated lipids, DHPC and CHAPS might form small bicelles at low q values that quickly grow in size as the relative amount of detergent is reduced. The presence of negatively charged lipids in the mixtures had an overall minor effect on the R H .
SAXS measurements were performed to gain further insight into the structural characteristics of the particles studied by DLS. As for DLS, SAXS measurements were taken at 25 • C and the R g values extracted from the SAXS curves are shown in Figure 2, lower panels. The R g data obtained from the model-independent analysis of the SAXS data were in close agreement with the DLS-based R H data, providing a similar picture of the q-dependent size changes in the particles.
The SAXS data were fitted to the elliptical bicelle model [17] implemented in SASview (http://www.sasview.org/, accessed on 1 October 2022) [19], consisting of a phospholipid bilayer modeled as an elliptical cylinder subdivided into three sections: heads, tails, and heads. The heads are symmetrical and are referred to as the disk face. The disk is surrounded by a belt or rim, representing the detergent's hydrophilic moiety. The hydrophobic part of the detergent is assumed to be incorporated into the tail region of the disk. The SAXS data were fitted by fixing the scale factor to the calculated volume fraction of the components in the solution obtained by summing the contribution of the molecular volumes of each component adjusted by its concentration. The electron densities of the components were initially set to the calculated values (Table S1) and allowed to vary to account for solvent penetration effects, while the thickness of the face, core, and rim were fitted together with the size of the cylinder's major radius and the ratio between the major and minor radii. The best fits to the data (solid lines) are included in Figures 3, 4A and S2, and the parameters derived from such fits are shown in Figures 4B, 5B and S7, and the values are reported in Table S2.
The bicelle model was found to produce adequate fits for all DMPC ( Figure 3) and DMPC/DMPG ( Figure S2) mixtures with the three detergents studied. The bilayer head groups (face) and core thickness values were comparable for bicelles composed of DHPC and CHAPS, and ranged around 23-25 Å for the tail region and 9-12 Å for the head groups, except for DMPC-CHAPS at q = 0.5. The headgroup and core thickness values are consistent with those of a DMPC lipid bilayer core in the fluid phase [20][21][22]. The rim thickness stabilized above q = 1.5 at approximately 8 and 4 Å for CHAPS and DHPC, respectively. The disks became more elongated with increasing q, growing by a factor of 1.6-2 from q = 0.5 to 5.
particle size and morphology for DMPC-detergent mixtures at the intermediate value of q = 2.5 was confirmed using cryo-EM ( Figure S3). Small spherical particles were observed for sodium cholate and slightly elongated particles were seen for CHAPS and DHPC, with the largest particles formed by DMPC-DHPC mixtures, in agreement with the SAXS and DLS measurements. The SAXS data and best fits for the mixtures composed of POPC and DHPC are shown in Figure 4. SAXS data for the other combinations of unsaturated phospholipids and detergents were also collected and are shown in the Supplementary Material (Figures S4-S6); however, these measurements did not follow a clear q-dependent trend and could not be adequately fitted using the bicelle model employed for the other lipid-detergent mixtures. Our SAXS analysis rationale was based on the q-dependent evolution of the scattering curves, which shifted gradually at lower Q and lower intensities as the detergent concentration was lowered. This suggested a continuous evolution of the bicellar structure, which was reflected in the progressive increase in the radius ( Figure 5). The absence of a clear trend in the unsaturated lipid data (except for the POPC-DHPC system) The data for the DMPC-sodium cholate mixtures clearly differed from those of the other DMPC-detergent mixtures. These particles were not only smaller in radius, but the bilayer core was significantly thinner compared with the particles formed with the other two detergents. The face and rim thickness were similar to those found for the DHPC and CHAPS mixtures. Thus, the overall shape of the DMPC-sodium cholate particle remained rather spherical, even at the highest q examined here, if we consider the total bilayer thickness (~10 × 2 +~19 =~49 Å) and the total disk diameter (~4 × 2 +~25 × 2 =~58 Å), together with the ellipticity value remaining close to 1. Therefore, it is plausible that DMPC-cholate mixtures form mixed micelles, rather than bicelles within this q range. The particle size and morphology for DMPC-detergent mixtures at the intermediate value of q = 2.5 was confirmed using cryo-EM ( Figure S3). Small spherical particles were observed for sodium cholate and slightly elongated particles were seen for CHAPS and DHPC, with the largest particles formed by DMPC-DHPC mixtures, in agreement with the SAXS and DLS measurements.
DHPC particles are shown in Figure 4 and Supporting Information Figure S7. The values obtained for the thickness of the head (face), core, and rim regions were comparable to those found for the corresponding saturated system. As indicated by the RH and Rg values, however, a significant increase in the size took place as q increased beyond 1. This was reflected in an increase in the major radius of the particles, which reached ~950 Å at q = 5.  The data for the DMPC-sodium cholate mixtures clearly differed from those of the other DMPC-detergent mixtures. These particles were not only smaller in radius, but the bilayer core was significantly thinner compared with the particles formed with the other two detergents. The face and rim thickness were similar to those found for the DHPC and CHAPS mixtures. Thus, the overall shape of the DMPC-sodium cholate particle remained rather spherical, even at the highest q examined here, if we consider the total bilayer thickness (~10 × 2 +~19 =~49 Å) and the total disk diameter (~4 × 2 +~25 × 2 =~58 Å), together with the ellipticity value remaining close to 1. Therefore, it is plausible that DMPC-cholate mixtures form mixed micelles, rather than bicelles within this q range. The particle size and morphology for DMPC-detergent mixtures at the intermediate value of q = 2.5 was confirmed using cryo-EM ( Figure S3). Small spherical particles were observed for sodium cholate and slightly elongated particles were seen for CHAPS and DHPC, with the largest particles formed by DMPC-DHPC mixtures, in agreement with the SAXS and DLS measurements. The ellipticity ratio defines the minor radius orthogonal to the major radius in the bilayer plane (B) Parameters of the best bicelle fits for the SAXS curves shown in Figure 3. Same colour scheme as in Figure 2. The full list of values is given in Table S2.
The increase in the parameter describing the major radius was accompanied by a concomitant decrease in the ellipticity ratio. Notably, the resulting minor radius remained between 20 and 30 Å for all q values, suggesting that the structures grew almost exclusively in one direction into highly elongated particles that resembled rods more closely than large, flat disks, as indicated by the relatively constant values of the minor radius and bilayer thickness across the q ratios explored. With the increasing size, a Bragg peak at a Q value of 0.096 Å −1 became apparent at q = 2.5, which increased in intensity at q = 3 and at q = 5 ( Figure  4C and Supporting Information Figure S8). Under the latter conditions, a faint second Bragg peak also started to appear at 0.192 Å −1 , which was twice the Q value of the first one. The primary Bragg peak corresponded to repeating structures with a correlation distance of ~65 Å. Given the elongated nature of the structures, it is reasonable to assume that this distance corresponded to the distance between the centers of two adjacent, side-by-side rods, resulting in a separation of 15-20 Å between their hydrophilic external regions. Bragg peaks did not appear in any other lipid detergent combination measured here and were also absent from the POPC-POPG DHPC sample, indicating that the presence of 20% POPG in the mixture was already sufficient to change the internal organization.
Finally, an attempt was made to produce bicelles using natural phospholipids from polar extracts of E. coli ( Figure 6). DLS data were collected for mixtures with DHPC, Figure 5. (A) Cartoon representing the bicelle model implemented in SASview used to fit the SAXS data. The ellipticity ratio defines the minor radius orthogonal to the major radius in the bilayer plane (B) Parameters of the best bicelle fits for the SAXS curves shown in Figure 3. Same colour scheme as in Figure 2. The full list of values is given in Table S2.
The SAXS data and best fits for the mixtures composed of POPC and DHPC are shown in Figure 4. SAXS data for the other combinations of unsaturated phospholipids and detergents were also collected and are shown in the Supplementary Material (Figures S4-S6); however, these measurements did not follow a clear q-dependent trend and could not be adequately fitted using the bicelle model employed for the other lipid-detergent mixtures. Our SAXS analysis rationale was based on the q-dependent evolution of the scattering curves, which shifted gradually at lower Q and lower intensities as the detergent concentration was lowered. This suggested a continuous evolution of the bicellar structure, which was reflected in the progressive increase in the radius ( Figure 5). The absence of a clear trend in the unsaturated lipid data (except for the POPC-DHPC system) indicated that the change in the structure was of a more complex nature, and likely arose from the coexistence of non-homogeneous structures affected by the lipid:detergent ratio in a non-trivial manner. Such heterogeneous structural composition is no longer represented by the simple single-core-shell particle model or bicelle model. The only combination that resulted in a clear trend was the POPC-DHPC series, which was, therefore, modeled with the bicelle form factor. The parameters for the best fits obtained from the POPC-DHPC particles are shown in Figure 4 and Supporting Information Figure S7. The values obtained for the thickness of the head (face), core, and rim regions were comparable to those found for the corresponding saturated system. As indicated by the R H and R g values, however, a significant increase in the size took place as q increased beyond 1. This was reflected in an increase in the major radius of the particles, which reached~950 Å at q = 5.
The increase in the parameter describing the major radius was accompanied by a concomitant decrease in the ellipticity ratio. Notably, the resulting minor radius remained between 20 and 30 Å for all q values, suggesting that the structures grew almost exclusively in one direction into highly elongated particles that resembled rods more closely than large, flat disks, as indicated by the relatively constant values of the minor radius and bilayer thickness across the q ratios explored. With the increasing size, a Bragg peak at a Q value of 0.096 Å −1 became apparent at q = 2.5, which increased in intensity at q = 3 and at q = 5 ( Figure 4C and Supporting Information Figure S8). Under the latter conditions, a faint second Bragg peak also started to appear at 0.192 Å −1 , which was twice the Q value of the first one. The primary Bragg peak corresponded to repeating structures with a correlation distance of~65 Å. Given the elongated nature of the structures, it is reasonable to assume that this distance corresponded to the distance between the centers of two adjacent, side-byside rods, resulting in a separation of 15-20 Å between their hydrophilic external regions. Bragg peaks did not appear in any other lipid detergent combination measured here and were also absent from the POPC-POPG DHPC sample, indicating that the presence of 20% POPG in the mixture was already sufficient to change the internal organization.
Finally, an attempt was made to produce bicelles using natural phospholipids from polar extracts of E. coli ( Figure 6). DLS data were collected for mixtures with DHPC, cholate, and CHAPS. Small particles were formed for DHPC at q = 0.5 (R H < 10 nm), while a broad particle size distribution was found for larger q values, ranging within 15 nm < R H < 100 nm. Mixed micelles and vesicles or worm-like micelles were most likely formed at low and high q values, respectively. Additional studies with detailed DLS measurements and EM image acquisition would be necessary to identify the exact structures formed, but this is beyond the scope of the present work. For mixtures with sodium cholate, on the other hand, small particles with R H < 10 nm were found for all q values. This suggests that small mixed micelles were formed in most cases, as was the case for the other lipids studied. Finally, a broad size distribution was found for mixtures with CHAPS ranging from a few nanometers in size (for q < 1.5) to approximately 100 nm for q > 2. This suggests that either small mixed micelles or vesicles/worm-like micelles were formed in this case as well.
Molecules 2022, 27, x FOR PEER REVIEW 8 cholate, and CHAPS. Small particles were formed for DHPC at q = 0.5 (RH < 10 nm), w a broad particle size distribution was found for larger q values, ranging within 15 n RH < 100 nm. Mixed micelles and vesicles or worm-like micelles were most likely form at low and high q values, respectively. Additional studies with detailed DLS meas ments and EM image acquisition would be necessary to identify the exact struct formed, but this is beyond the scope of the present work. For mixtures with sodium late, on the other hand, small particles with RH < 10 nm were found for all q values. suggests that small mixed micelles were formed in most cases, as was the case for the o lipids studied. Finally, a broad size distribution was found for mixtures with CH ranging from a few nanometers in size (for q < 1.5) to approximately 100 nm for q > 2.
suggests that either small mixed micelles or vesicles/worm-like micelles were forme this case as well.

Discussion
The saturated phospholipid mixtures with DHPC allowed bicelles in the q ra studied to be obtained at 25 °C in TBS, except for q = 0.5 (Figures 2-5). Bicelles were obtained for saturated phospholipid mixtures with CHAPS. This was established by SAXS measurements (Figure 5), which demonstrated a clear bilayer structure with s metric heads and a lipid core having layer thicknesses consistent with the values expe for fluid, DMPC bilayers [20,22] when using saturated phospholipid mixtures w CHAPS or DHPC. For q = 0.5, the core thickness was lower than 20 Å and the disk r were 30 and 20 Å for DHPC and CHAPS, respectively. This suggests that mixed mice rather than disks, were formed in this case. This is in agreement with a previous stud DMPC-DHPC mixtures at q < 1 [18]. CHAPS is a detergent often used in membrane tein purification [23] as it maintains the functionality of the protein of interest. Ea NMR studies also found strong indications that bicelles could successfully be formed

Discussion
The saturated phospholipid mixtures with DHPC allowed bicelles in the q range studied to be obtained at 25 • C in TBS, except for q = 0.5 (Figures 2-5). Bicelles were also obtained for saturated phospholipid mixtures with CHAPS. This was established by the SAXS measurements ( Figure 5), which demonstrated a clear bilayer structure with symmetric heads and a lipid core having layer thicknesses consistent with the values expected for fluid, DMPC bilayers [20,22] when using saturated phospholipid mixtures with CHAPS or DHPC. For q = 0.5, the core thickness was lower than 20 Å and the disk radii were 30 and 20 Å for DHPC and CHAPS, respectively. This suggests that mixed micelles, rather than disks, were formed in this case. This is in agreement with a previous study of DMPC-DHPC mixtures at q < 1 [18]. CHAPS is a detergent often used in membrane protein purification [23] as it maintains the functionality of the protein of interest. Earlier NMR studies also found strong indications that bicelles could successfully be formed by using DMPC and CHAPS [14]. Therefore, CHAPS can, with confidence, be used instead of DHPC to form bicelles containing membrane proteins in the q-range of 1 < q < 5. Indeed, it was previously found that the use of CHAPS increased the stability of the transmembrane protein Opsin in DMPC-detergent bicelles [24]: this protein is particularly unstable in detergent or in DHPC-based bicelles.
Sodium cholate is another detergent often used in membrane purification. Mixtures of saturated phospholipids with sodium cholate exhibited no significant dependency of the lipid-detergent q range on particle size ( Figure 2). Moreover, the lipid core was thinner than expected for a fluid DMPC bilayer ( Figure 5). This suggests that mixed micelles are formed, rather than bicelles, when using sodium cholate as a detergent. Accordingly, cryo-EM corroborated the lack of a disk-like appearance for DMPC-cholate mixtures as compared with DMPC-DHPC or CHAPS mixtures at q = 2.5 ( Figure S3).
For unsaturated phospholipids, bicelles were only found under very limited experimental conditions (for q = 0.5 and 1). Higher q values produced particles with significantly larger sizes. This is expected, as DMPC-DHPC bicelles tend to merge and grow into vesicles at temperatures well above the gel-fluid phase-transition temperature, T m , for DMPC [25] or upon further increasing q [26]. The T m of POPC and POPG is around 0 • C. Therefore, these lipids are in the fluid phase at room temperature, and it is expected that elongated disks or vesicles are formed at intermediate q values for both DHPC and CHAPS. Indeed, such structures have been observed previously for DMPC-DHPC at q = 3-4 measured in 20 mM HEPES buffer at pH 7.4 enriched with 150 mM NaCl [27]. Detailed microscopy studies would enable whether elongated bicelles or vesicles are formed at intermediate and high q values for unsaturated phospholipids to be determined. Further determination of the non-bicellar structures formed by the mixtures hereby presented beyond the scope of the present work.
Polar E. coli extracts have previously been shown to produce bicelles at q = 0.5 when mixed with DHPC, as studied by NMR, yielding R H = 39 ± 5 Å [28]. The authors suggested that the larger size (and variation) found in R H was due to the longer lipid tails present in the E. coli bicelles as compared with the traditional DMPC/DHPC bicelles. Indeed, the authors found that the composition of the bicelles depended on whether the lipid extracts were taken from cells grown at 20 or 37 • C, with major differences found in the main lipid types present, as well as the acyl length. This could be related to the natural extract mixture being above or below their respective T m , although this transition is rarely sharp in natural lipid extracts [29]. For POPC/POPG-DHPC at q < 1, R H is~40 Å (Figure 2), while, for the polar E. coli extracts, R H is~8 and 50 Å for q = 0.5 and 1, respectively ( Figure 6). Thus, our data suggest that bicelles produced with polar E. coli lipid extracts could potentially only be obtained at q = 1 for DHPC. However, the polar E. coli extract-DHPC mixtures exhibited broad size distributions with no clear trend with q, as measured by DLS ( Figure 6); thus, it is not certain whether the structures observed correspond to bicelles at all. The use of cholate or CHAPS as detergents with polar E. coli lipid extracts did not show any promising results regarding the formation of bicelles.

Bicelles Preparation
For the preparation of the lipids and detergent stock solutions in organic solvents, phospholipids and DHPC were dissolved in chloroform, sodium cholate in methanol, and CHAPS in a 4:1 methanol:chloroform mixture. The latter required sonication and heating to 40 • C to completely dissolve. Lipid films were prepared by mixing phospholipids and detergents at the desired q to obtain a final phospholipid concentration of 5 mg/mL in all samples. Lipids and detergents were mixed in organic solvents at seven different lipid:detergent molar ratios to obtain q values of 0.5, 1, 1.5, 2, 2.5, 3, and 5, and the samples were dried overnight using a miVac Quattro concentrator (Genevac Ltd., Ipswich, UK). The dry lipid:detergent films were resuspended at a lipid concentration of 5 mg/mL in TBS buffer at pH 7.4 containing the desired detergent at a concentration just above the nominal critical micelle concentration (CMC). The detergent concentrations in the buffers used were 5 mM, 10 mM, and 16 mM for CHAPS, DHPC, and sodium cholate, respectively. The samples were subjected to three cycles of freeze-thawing using liquid nitrogen and a water bath heated to 45 • C, and vortexed for 1 min after each cycle. Samples were then centrifuged at 13,000 rpm for 10 min and the supernatants were used for further analysis.

Dynamic Light Scattering (DLS)
The size of the different bicelle preparations was determined using multi-angle dynamic light scattering (MADLS) on a Zetasizer Ultra (Malvern Panalytical, Ltd., Malvern, UK) using a quartz cuvette (ZEN2112) with a final volume of 90 µL. All measurements were performed on freshly prepared bicelle samples, mostly in duplicate or triplicate, at 25 • C and data were processed using the software Zetasizer Ultra-Pro ZS Xplorer v1.31.

Small-Angle X-ray Scattering (SAXS)
SAXS measurements were performed at beamline P12 operated by EMBL Hamburg at the PETRA III storage ring (DESY, Hamburg, Germany). Measurements were performed using an X-ray beam of 10 kEV, corresponding to a wavelength of 1.24 Å, a beam size of 250 × 100 µm 2 (FWHM), and a flux of 5.10 12 photons per second.
Samples and buffers were loaded in a flow-through capillary using the beamline sample changer and flown during X-ray exposure to limit radiation damage. Scattered photons were collected on a Pilatus 6M detector (Dectris, Baden, Switzerland). With a sample-to-detector distance of 3 m, the SAXS intensities can be collected for Q ranging between 0.0027 Å −1 and 0.74 Å −1 .
Data were also collected on the buffer solution used to resuspend the bicelles. These data contained the scattering of the buffer and of the instrumental background that can be subtracted from the sample curve to obtain the bicelle form factor. Initial data reduction was performed using the data analysis pipeline SASFLOW [30]. For each measurement, 40 successive frames with exposure of 100 ms were collected. The images obtained were radially average. The successive frames were compared, and those affected by radiation damage were discarded. The remaining frames were averaged, and the corresponding buffer curve was subtracted to obtain the form factor of the micelles.
Data obtained using the laboratory SAXS instrument (XEUSS3.0, Xenocs, Grenoble, France) at Malmö University showed consistent results on replicate samples, albeit with a slightly reduced Q range and with reduced counting statistics. The agreement between the synchrotron and laboratory SAXS supports the stability and reproducibility of the bicelle preparation, and the absence of any effects of radiation damage.

SAXS Analysis
Calculations of the radius of gyration (R g ) were performed using the AUTORG function [31] available in the Primus software version 3.0, EMBL, Germany [32].
SAXS data were fitted using SasView 4.2.2 [19]. Bicelles were modeled with the form factor of a core-shell elliptical cylinder described in the software documentation [33]. Briefly, the hydrophobic core of the bicelles was modeled as a cylinder with an elliptical base defined by a major and minor radius. The hydrophilic lipid head groups were described by two symmetric cylinders adjacent to the bottom and top faces of the central cylinder, while the detergent ring surrounded it, covering the length of the hydrophobic core, as described in Figure 5. The scale factor that defined the intensity of the simulated data was set to match the volume fraction of scattering components calculated by adding together the molecular volumes of the species in solution according to their concentrations (Table S1) and allowed to vary within a 10% range from its calculated value. The same molecular volumes were used to estimate the electron densities of the components to define the starting values of the fits. The electron density of the hydrophobic core was kept fixed to its calculated value, while that of the head group was allowed to vary to consider solvent penetration within the hydrophilic regions.

Cryogenic Electron Microscopy (Cryo-EM)
Cryo-EM image acquisition was performed on a JEM-2200FS transmission electron microscope (JEOL, Tokyo, Japan) at the National Center for High-Resolution Electron Microscopy (nCHREM) at Lund University. The microscope was equipped with a fieldemission electron source, a cryo pole piece in the objective lens, and an in-column energy filter (omega filter). Zero-loss images were recorded at an acceleration voltage of 200 kV on a TemCam-F416 camera (TVIPS) using SerialEM under low-dose conditions. Samples were prepared using an automatic plunge freezer system (Leica EM GP, Wetzlar, Germany) with the environmental chamber set to 21 • C and 90% relative humidity. A 4 µL droplet of the sample solution was deposited on a lacey formvar carbon-coated grid (Ted Pella, Redding, CA, USA) and was blotted with filter paper to remove excess fluid. The grid was then plunged into liquid ethane (around −183 • C) to ensure rapid vitrification of the sample in its native state. The specimens were then stored in liquid nitrogen (−196 • C) and, prior to imaging, transferred to the microscope using a cryo transfer tomography holder (Fischione, Export, PA, USA, Model 2550).

Conclusions
This study clearly shows that bicelles can be obtained with saturated phospholipids (DMPC and DMPG) using both DHPC and CHAPS as detergents for q > 1. Mixed micelles were likely formed for q < 1. The study further demonstrates that no bicelles can be obtained when using sodium cholate as a detergent, regardless of the q used, and instead very small mixed micelles are formed. In the case of unsaturated phospholipids (POPC and POPG), bicelles can be obtained by using both DHPC and CHAPS as detergents, but only for q < 1 Finally, our results show that bicelles cannot be successfully produced with E. coli polar extracts, at least with the detergents used in this study and within the q range used. Note that the experimental conditions in this study were saline TBS buffer at pH 7.4 and 25 • C, which are physiologically relevant conditions suitable for the study of transmembrane and membrane-bound proteins.
Supplementary Materials: The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/molecules27217628/s1, Table S1: Chemical composition, molecular volumes, and electron densities of the lipids and detergent components. a is from [34], b is from [35], c is from [36], and d is estimated using the VegaZZ software [37]. Table S2. Complete list of the values describing the fits obtained for the curves shown in Figure 3, Figure 5 and Figure S2 and plotted in Figure 4B, Figure 5B and Figure S7. Figure S1: DLS size distributions of the lipid-detergent mixtures measured in this study; Figure S2: SAXS curves and best fits to the data for the systems containing DMPC-DMPG (80:20) and the three detergents. Figure S3: Cryo EM images of DMPC mixed with cholate (left), CHAPS (center), and DHPC (left) at a q ratio of 2.5; scale bar = 50 nm. The schematic drawings underneath represent the structures obtained from the fitting of the SAXS data drawn to scale using the thickness parameters obtained from the fits. Figure S4: SAXS curves for the POPC-POPG (80:20) samples mixed with DHPC. Figure


Unravelling the landscape of skin cancer through single-cell transcriptomics

Highlights • Human skin is a multifunctional and heterogeneous organ with a vast cellular pool.• Single-cell RNA-sequencing provides an unprecedented opportunity to investigate and understand the heterogeneous landscape of skin cancers.• Analysis of single-cell skin cancer transcriptomes opens new avenues of investigation, including the potential for precision-targeted therapy.

Introduction
Cellular diversity and heterogeneity underlie the skin's many functions, which include sensation, stabilization of body temperature, maintaining fluid and electrolyte balance, vitamin and hormone synthesis, protection from external insults, and immunologic defense [1][2][3][4]. The skin is divided into three broad compartments-the epidermis, the dermis, and hypodermis/subcutaneous layer. Each compartment contains specific anatomical structures and diverse cell populations that contribute to skin physiology and homeostasis [1][2][3][4][5][6]. The epidermis is arranged in stratified epithelial layers that establish the skin permeability barrier. While keratinocytes comprise >90% of the epidermal cell population, the epidermis also contains melanocytes, Langerhans cells, and Merkel cells [1][2][3][4][5][6]. By comparison, the dermis is normally paucicellular and characterized by an extracellular matrix of collagen and elastin fibers that provides structural support as well as a niche for sparsely distributed fibroblasts and a variety of immune cells [1][2][3][4]7]. Skin adnexa, blood and lymphatic vessels, and nerve bundles are also housed in the dermis [1][2][3][4]. The hypodermis is primarily composed of adipocytes organized into fat lobules and loose connective tissue; however, it also contains fibroblasts and macrophages (Fig. 1). Crosstalk between different cell populations, such as structural epithelial cells and immune cells, is increasingly appreciated to play an important role in both epidermal homeostasis as well as disease states [1,3,4].
Our aging population, rising global temperatures, and challenges to changing sun-related attitudes and behaviors all contribute to the rising incidence of skin cancer. Keratinocyte carcinomas, namely, basal cell carcinomas (BCC) and cutaneous squamous cell carcinomas (cSCC) [8], account for the majority of these cases and are the most common cancers in the United States [9,10]. While less common, melanoma and Merkel cell carcinoma (MCC) also demonstrate growing incidence rates. Histopathologically, these tumors are characterized by nests or sheets of neoplastic cells within the epidermis and/or involving the dermis as well immune cell infiltrates that reflect the host immune response. Local treatment of small and/or early skin cancers by surgical excision or radiation therapy is often curative; however, effective treatments for advanced or metastatic tumors have met with more measured success. The heterogeneity of tumor cells and cells comprising the tumor microenvironment (TME) often underlies therapeutic resistance, suggesting that a deeper understanding of these cell populations and their interactions might facilitate the development of more effective targeted treatments. Single-cell RNA-sequencing (scRNA-seq) now provides an unprecedented opportunity to study malignant cell heterogeneity and TME complexity at single-cell resolution. A number of recent review articles details how scRNA-seq has recently been applied to study skin homeostasis as well as cutaneous inflammatory diseases and cancers [8,11,12]. Here we discuss how our understanding of skin cancer biology has advanced using bulk transcriptomics as well as single-cell approaches and highlight different methods of scRNA-seq analysis.

Recent applications of RNA-seq in skin cancer
When compared to hybridization-based microarrays, nextgeneration RNA sequencing (RNA-seq) possesses several strengths, including greater dynamic range and the capacity for de novo transcriptome assembly. In the skin cancer setting, whole transcriptome profiling by RNA-seq has been used to identify novel genes, such as noncoding RNAs and lowly expressed mRNA isoforms/transcripts. More recently, RNA-seq-based approaches have also been used to predict treatment response, nominate new therapeutic targets, and advance our understanding of disease progression as well as pathogenesis. For example, a recent analysis of melanocytic nevi and primary melanomas by Kunz et al. using RNA-seq revealed separate transcriptome signatures, with enrichment of genes driving BRAF/MEK inhibitor resistance in NRAS-mutant nevi and primary melanomas, while genes associated with resistance to PD-1 inhibition were expressed in lesions with wild type NRAS [13]. This work defined two different trajectories for melanoma development, suggesting that treatment resistance is determined early in melanomagenesis [13]. Additional work by Svedman et al. using targeted RNA sequencing (Ion AmpliSeq) to characterize melanomas following immune checkpoint inhibition similarly demonstrated enrichment of key pathways such as DNA replication, chromatin remodeling, and cell cycle, which may predict long-term response to therapy [14]. RNA-seq has also been used to highlight novel targets and signaling pathways in BCC. Wan et al. analyzed the transcriptomes of BCC with high tumor purity compared to patient-matched normal skin and identified enrichment of C2H2-type Zinc Finger genes, including Gli transcription factors as well as many others not functionally characterized in this malignancy [15]. This study, along with a larger pooled analysis of publicly available RNA-seq data by Litvinov et al. confirmed known dysregulated pathways in BCC and also suggested new therapeutic targets in the Wnt/β-catenin as well as IL-17 signaling pathways [15,16]. Other groups have used RNA-seq to elucidate the molecular underpinnings of skin malignancies. Using human specimens that represent increasingly dysplastic states on the cSCC spectrum as well as the equivalent tissues in SKH-1E cSCC-prone mice, Chitsazzadeh et al. distilled four TFs -ETS2, SP1, FOXF2 and AP1 -that appear to regulate the continuum of cSCC development [17]. Similar results were obtained by Das Mahapatra et al., who compared cSCC to unmatched normal skin controls and also reported enrichment of genes regulating immunological pathways, highlighting the importance of immunotherapy in this malignancy [18]. In MCC, which is divided into virus-positive and virus-negative categories based on the presence or absence of Merkel cell polyomavirus (McPyV) transcripts, Starrett et al. used RNA-seq to compare virus-positive to virus-negative MCC transcriptomes and demonstrated that genes regulating the cell cycle were increased in the former while DNA repair genes were downregulated in the latter [19]. These findings demonstrated a role for McPyV in controlling the MCC transcriptome and provide an explanation for the high mutation burden observed in virus-negative tumors [19].
As demonstrated by the aforementioned studies, bulk RNA-seq-based analyses are a valuable mainstay of investigative skin cancer biology. There are many advantages to using this approach, such as detection of mRNA isoforms and RNA species, high read capture capacity which enables detection of low abundance genes, low signal-to-noise ratio, and cost effectiveness. However, when applied to cancer, a few important limitations should be considered. First, details regarding the cellular context of the observed mRNA changes are not readily provided by bulk RNA-seq, although methods now exist to deconvolute this data into known cell populations using reference signatures. Second, bulk sequencing measures the average expression of a pooled population of cells, raising the possibility that transcriptomic changes of rare populations may be masked. This is true of many tissues including skin, where heterogeneity is present not only at the level of diverse cell types and states (Fig. 1), but also evolves from disease progression and treatment. The emergence of therapeutic resistance following cancer chemotherapy, for example, is a setting in which bulk sequencing can fall short due to its inability to distinguish between drug-sensitive and resistant clones. Additionally, the ability to explore RNA trajectory, cellspecific multi-omics or cell-cell communication is limited with bulk RNA-seq. Thus, these and similar questions are increasingly being addressed by single-cell approaches that provide a more granular view of cell types and states, potentially revealing new targetable transcriptomic changes.
The two widely used methods for single cell sequencing are 10X genomics and Smart-seq2. Smart-seq2 has been shown to detect low abundance transcripts and isoforms, while 10X genomics is better at detecting rare cell types [20]. Using these methods and scRNA-seq data analysis approaches, researchers have highlighted cellular/spatial context of skin cancer transcriptome and deciphered novel cell-to-cell communications.

Cellular context
A key goal of scRNA-seq is to characterize the cell types present within a sample, as this allows researchers to understand cell-type specific gene expression. An important first step when analyzing scRNA-seq data is to consider batch effects, or gene expression patterns caused by non-biological factors that can cause erroneous interpretation of cell types. Removing this technical noise for proper biological interpretation is referred to as integration. Methods like Seurat CCA, Harmony, scVI, scanorama, and scMC have been developed to integrate cells from multiple technical conditions [21][22][23][24][25]. Although integration methods help remove technical noise, when applied incorrectly they can overcorrect and remove biological variations.
Many approaches have been developed to identify cell types. Manual annotation is a popular technique combining dimensionality reduction and unsupervised clustering followed by manual inspection of marker genes. This approach works best for smaller studies or when most of the cells being analyzed are expected to have few known marker genes. For larger studies of samples with previously profiled cell types, automated classification methods can greatly speed up annotation. Some of these approaches, such as SingleR, CellAssign, and Garnett rely on transcriptomic reference profiles or gene signatures to label cell types [26][27][28][29][30]. By contrast other methods, including Azimuth, ProjecTILs, and scArches projects new cells to a reference atlas of previously annotated cell types [27,31,32]. Large consortium efforts such as Human Cell Atlas are currently cataloging major cell types and generating reference atlases that can be used to map and classify data from new studies [33].
After identifying broad cell types, further clustering is often performed on cells within a specific lineage (e.g., tumor cells, fibroblasts, Tcells) to identify subpopulations of cells with transcriptionally distinct states. This analysis is usually exploratory in nature and thus performed manually unless prior studies have already determined marker genes for distinct cell states. In the cancer setting, identification of cell types and subclusters is often used to characterize intra-and intertumoral heterogeneity and can also be applied to immune and stromal cells to study the transcriptional heterogeneity of the TME.
For example, in the first scRNA-seq study of skin cancer, Tirosh et al. analyzed 4645 cells isolated from 19 patients with melanoma and showed that across all tumors, malignant cells from the same melanoma contained two distinct transcriptional states not distinguishable by bulk RNA-seq that were characterized by either high expression of MITF or AXL. Cells with high AXL expression exhibited drug resistance to RAF/ MEK inhibition, suggesting they drive tumor recurrence [34]. This work also highlighted differential T-cell activation, expansion, and clonal exhaustion programs across different patients [34]. By sequencing single-cells from short-term culture of three melanomas with differing BRAF/NRAS genotypes, Gerber et al. derived gene expression modules from different cell subpopulations that included those with high MITF or AXL and also revealed targetable upregulation of the CDK4 and CDK2 cell cycle-dependent kinases responsive to palbociclib [35]. Yao and colleagues, identified similar patterns of intratumoral heterogeneity in BCC by utilizing scRNA-seq to characterize tumor cells based on expression of nuclear myocardin-related transcription factor (nMRTF), which correlates with Smoothened (SMO) inhibitor resistance and is targetable using AP-1 inhibitors [36]. In cSCC, Ji et al. also identified four different keratinocyte subpopulations, including tumor-specific keratinocytes (TSKs). TSKs are observed at the leading edge of the tumor and are thus positioned to act as a hub in cSCC tumor-stromal interactions [37]. Tumor subpopulation clustering was also used by Paulson et al. to elucidate mechanisms of treatment relapse in McPyV+ MCC by demonstrating tumor cells adapt to the immunological pressure created by autologous McPyV-specific CD8 + T cells and immune checkpoint inhibitors by suppressing the HLA specific to the targeted McPyV epitope [38].
When applied to cancer-associated immune and stromal cells, the identification of cell types and subclusters has broadened our understanding of the transcriptional heterogeneity of the TME. Several studies have examined the landscape of T-cells in skin cancers, finding a diverse range of cells including cytotoxic, effector, exhausted, regulatory, and helper T-cells. Yost et al. coupled scRNA and T-cell receptor sequencing (scTCR-seq) to study advanced BCC following PD-1 blockade and showed that novel T-cell clones infiltrate the TME of advanced BCC upon treatment, replacing their pre-existing exhausted counterparts rather than reinvigorating them [39]. Deng et al. performed a meta-analysis of 59 melanomas and catalogued the transcriptional states of CD8 + T-cells. They defined seven cytotoxic, exhausted, and naive/memory subpopulations, including a subset of exhausted T-cells associated with poor prognosis and characterized by high expression of PMEL, TYRP1, and EDNRB [40]. Frazzette et al. recently used scTCR-seq and gene expression data to compare the T-cell landscape in cSCC from skin cancer-prone immunocompromised organ transplant recipients (OTR) to that of immunocompetent individuals [41]. This effort highlighted several differences between the tumor-infiltrating lymphocytes (TILs) of these two groups, with a reduction in cytotoxic T-cell number, TCR clonotypes, and clonal expansion observed in OTR [41]. Other studies, such as that by Davidson and colleagues, have focused on characterizing cancer-associated fibroblasts (CAFs) [42]. Using the murine B16 melanoma model, these authors defined three functionally distinct CAF subpopulations that are positioned to play roles in immune crosstalk, extracellular matrix remodeling, and cytoskeletal reorganization [42].

Cell-cell communication
The availability of expression profiles from multiple cell types within a single sample allows detection of intercellular signaling via ligandreceptor expression analysis [43,44]. Several computational methods have been developed, including CellPhoneDB and CellChat [45,46]. These tools utilize databases with prior knowledge of ligand-receptor pairs or complexes to infer cell-cell communication based on expression of ligands and receptors from annotated cell groups (Fig. 2). This approach has proved to be especially useful for investigating interactions between tumor cells and their microenvironment as well as immune-stromal crosstalk in skin cancers (Fig. 2).
Within the TME, CAFs are well known to modulate tumor growth and progression. In their work identifying the TSK keratinocyte subcluster in cSCC, Ji and colleagues used NicheNet [37] to highlight TSK signaling to CAFs through ligand-receptor pairing of MMP9-LRP1 and TNC-SDC1 37 . Similarly, work by Davidson et al. used CellPhoneDB to identify a subset of CAFs in the murine melanoma TME that express the immunomodulatory C3 protein and are in close contact with C3aR-positive myeloid cells; disrupting this interaction reduced tumor growth and resulted in fewer myeloid cells that suppress CD8 + T-cells [42]. Deng et al. also applied CellPhoneDB to their re-analysis of publicly available single-cell melanoma transcriptomes to show that ligand-receptor interactions involving exhausted T-cells, such as CCL5-CCR5 and CD74-MIF, are central to communication between different CD8 + subpopulations [40]. Guerrero-Juarez et al. used CellChat to investigate tumor-stromal interactions in BCC and identified that the nearby fibroblasts act as an inflammatory hub regulating BCC growth through induction of heat shock proteins [47]. Knowledge of ligand-receptor pairing was also essential to work performed by Miao and colleagues, who performed adoptive cytotoxic T-cell transfer (ACT)-based immunotherapy in a mouse model of cSCC to investigate mechanisms of immune evasion [48]. Using scRNA-seq and lineage tracing, these authors identified a population of CD80 + tumor-initiating stem cells that are refractory to ACT and engage with CTLA4-expressing cytotoxic T-cells, reducing their activity and promoting tumor relapse [48].

Copy number alterations
Genetic variation is traditionally assayed using exome or genome sequencing, which provides data on mutations, structural variants, and genomic copy number. When applied to cancer, these methodologies allow researchers to query the clonal structure of tumors, estimate tumor purity and ploidy, and distinguish neoplastic from healthy samples. Single cell DNA-sequencing (scDNA-seq) enables elucidation of genetic heterogeneity at higher resolution than bulk exome/genome approaches; however, it is not used as widely as scRNA-seq due to its high cost, lack of commercial kits, and data quality issues [49]. Method developers have attempted to circumnavigate these issues by estimating single-cell copy number profiles from scRNA-seq data with tools such as inferCNV, honeyBadger, copyKat, and CaSpER [50][51][52]. These approaches have been used to distinguish neoplastic from normal cells and understand clonal lineages within tumors. For example, in their scRNA-seq analysis of melanoma, Tirosh et al. classified melanocyte-like cells as malignant, non-malignant, or intermediate based on a copy number alteration score [34]. Only malignant cells were included in subsequent analyses of tumor gene expression programs, thus avoiding contamination from normal melanocytes that can occur when analyzing bulk RNA-seq data [34]. Yost and colleagues demonstrated concordance between copy number variants inferred from single-cell transcriptomes and those called using whole exome sequencing data in BCC [39].

Gene trajectory
Many biological processes such as differentiation create cells that exist along a continuum rather than distinct transcriptional states. Trajectory inference methods have been developed to infer the relationship between cells across a continuous process [53,54]. RNA velocity analysis is a related technique that uses the proportion of unspliced and spliced mRNA molecules within a cell to provide temporal ordering along the trajectory [55] (Fig. 3). When applied to skin cancer, trajectory analysis has also been used to study tumor evolution. Wouters et al. employed Monocle-2 to characterize phenotype switching in melanoma [56]. Their findings demonstrate that melanoma cells transition from a melanocytic transcriptional state to a mesenchymal-like one through a stable intermediate state characterized by unique chromatin features and transcriptionally regulated by SOX6, NFATC2, EGR3, ELF1 and ETV4 56 . Another trajectory study by Su and colleagues found that BRAF V600E melanomas exhibit a bifurcated developmental trajectory in response to BRAF inhibition that is dependent on MITF expression [57]. Trajectory analyses have also provided new insight into the TME response to treatment. Using scRNA and T-cell receptor sequencing, Yost et al. showed that novel T-cell clones infiltrate the TME of advanced BCC following PD-1 blockade, replacing their pre-existing exhausted counterparts. The authors used trajectory analysis to analyze the evolution of T-cells in response to anti-PD-1 therapy and discovered they follow a bifurcated trajectory, becoming either terminally activated or exhausted, with the latter state marked by increased expression of PDCD1 and HAVCR2 39 . Deng and colleagues observed similar findings in melanoma, with their trajectory analysis suggesting that CD8 + T-cells convert to either exhausted or cytotoxic terminal states [40].

Bulk deconvolution
While scRNA-seq offers new opportunities to study skin cancer, bulk RNA-seq remains more common due to its lower cost and ease of use. Transferring findings from single-cell experiments to bulk datasets to take full advantage of available data resources is thus an important part of single-cell analysis. The existence of several publicly available skin cancer datasets with bulk microarray or RNA-seq data allows researchers to validate findings from single-cell experiments in larger cohorts of bulk data such as The Cancer Genome Atlas (TCGA) [58]. One common approach is to derive cell-type or tumor subpopulation-specific gene signatures from single-cell data and then score bulk data using these signatures, allowing signature activity to be correlated with clinical outcomes such as overall survival or metastasis [34]. More advanced deconvolution algorithms such as CIBERSORTx, MuSiC, and others use single-cell or bulk reference profiles to estimate cell-type abundance and gene expression from bulk tumor samples [59][60][61]. In the skin cancer setting, deconvolution has successfully been applied to study tumor-stromal interactions. Tirosh and colleagues used gene signatures derived from single-cell melanoma transcriptomes to deconvolute 471 TCGA melanoma samples and showed that CAF abundance correlates with AXL expression [34]. They also identified a set of CAF-expressed genes that correlates with T-cell infiltration [34]. Similarly, after deriving signatures for T-cell exhaustion and cytotoxicity from single-cell data and scoring TCGA melanoma samples with these signatures, Deng et al. identified cytotoxic and exhausted subpopulations that predict prognosis [40]. Other studies have used deconvolution algorithms to link levels of infiltrating immune cells with survival and nominate genes that may be involved in this process [62][63][64]. Deconvolution was also been applied to cSCC by Ji et al., who linked TSK-like expression with CAF activity in several TCGA cancer types [37].

Spatial context
Single-cell approaches are increasingly being combined with spatial transcriptomics (ST) to gain further insight into the complex nuances of different cell types and states. By providing tissue context, ST enables researchers to understand gene expression as it relates to cell position within a tissue. The two widely used methods for ST are imaging-based spatially resolved transcriptomics and in situ barcoding-based spatially resolved transcriptomics. Imaging-based spatially resolved transcriptomics methods such as MERFISH, seqFISH and ABER-FISH relies on multiplexed fluorescence in situ hybridization (FISH) and expansion microscopy to detect several RNA molecules [65][66][67][68]. In situ barcoding-based spatially resolved transcriptomics methods such as 10X Visium, Slide-seq and high-definition spatial transcriptomics (HDST) utilizes DNA barcoding to map out gene expression [68][69][70][71]. Although these methods can interrogate the entire transcriptome, they are limited by their tissue resolution, which varies from 55 μm to 2 μm. ST methods have been applied to skin cancer to localize tumor subpopulations and may shed light on their role in cancer progression (Fig. 4). Using a highly multiplexed and sensitive platform called CEL-Seq2 Baron and colleagues used scRNA-seq to characterize zebrafish melanomas and identified a tumor subpopulation that expresses a stresslike transcriptional program [72]. ST was used to demonstrate that these cells are enriched in tumor regions, but not in the normal surrounding tissue [72]. After confirming the presence of stress-like cells in human melanomas and other cancer types, the authors showed that the stressed state efficiently seeds new tumors in zebrafish and is associated with resistance to MEK and BRAF inhibitors. In another study, Thrane et al. used the 10X Visium ST platform to characterize melanoma lymph node metastases, highlighting differential gene expression profiles as well as multiple melanoma signatures across the same region in a single tumor biopsy [73]. Poor patient survival was associated with greater expression heterogeneity in the area of transition between melanoma and lymphoid tissue [73]. Ji et al. similarly used the 10X Visium ST platform to demonstrate that cSCC TSKs and basal cells localize to the tumor leading edge and confirm the presence of vascular and CAF-enriched transcripts in the surrounding fibrovascular niche [37]. By incorporating a novel multi-parameter tissue imaging workflow named Pick-Seq, Nirmal et al. showed that invasive cutaneous melanoma possesses a unique cellular microenvironment at the tumor-stromal boundary that harbors different cell types supported by cytokine gradients [74].
Some of the limitations of ST are that it lacks single-cell resolution and the read quality as well as the number of transcriptomes captured is often tissue-dependent. HDST notably enables users to reach a resolution of 2µM [68,70], providing the possibility of subcellular resolution. This method currently requires specialized analytics as well as bioinformatic expertise, and is limited by low sensitivity of mRNA capture [75]. Moreover, both scRNA-seq and ST do not consistently detect genes with lower expression. To overcome these barriers, approaches such as RNA in situ hybridization, which offers a targeted but more sensitive alternative, can be used to generate spatial information at single-cell resolution. A recent study by Tran et al. developed a multimodal strategy, Spatial Transcriptomic and RNA in situ Hybridization (STRISH), to study cancer-immune cell crosstalk at the genome-wide level in keratinocyte cancers [76]. This approach starts by inferring ligand-receptor interactions in the cancer-immune cell context using scRNA-seq and ST. Cell-cell interactions are then visualized by RNA in situ hybridization and quantitated using digital droplet PCR [76]. Using this analysis pipeline, the authors demonstrated co-expression of the IL34-CSF1R and THY1-ITGAM ligand-receptor pairs in cancer cell nests and areas of immune cell infiltrate [76]. We predict that similar strategies combining different technologies to overcome the limitations of each approach in isolation will be increasingly adopted to discover and confirm critical cell-cell interactions in the cancer setting.

Conclusions
Bulk RNA-seq remains a valuable method for defining the coding and non-coding transcriptome; however, when applied to the skin and its attendant malignancies, scRNA-seq has revealed robust expression differences between cell types and states, providing new insight into the cellular diversity and heterogeneity that occurs in these settings (Table 1). While this approach is accompanied by its own limitations, most notably dropout events and the scalability of existing data analysis methodologies, the aforementioned studies demonstrate how scRNA-seq has enabled an improved understanding of molecular events regulating skin cancer progression. We predict that future studies will utilize singlecell transcriptomics to interrogate larger cohorts of treatment-resistant and relapsing skin cancers, with an eye towards predicting therapy response. As sequencing costs fall, scRNA-seq may also increasingly be used in the clinical setting to devise targeted personalized treatment strategies for advanced skin cancers that cannot be managed surgically.

Declaration of Competing Interest
None.


s2fieldsofstudy: Biology Medicine

text:
Exploration the global single-cell ecological landscape of adenomyosis-related cell clusters by single-cell RNA sequencing

Background: Adenomyosis (AM) is a common benign uterine disease that threatens the normal life of patients. Cells associated with microenvironmental immune ecology are crucial in AM, although they are not as well understood at the cellular level. Methods: Single-cell sequencing (scRNA-seq) data were used to construct an AM global single-cell map, to further identify relevant cell clusters and infer chromosomal copy number variation (CNV) in AM samples. The biological functions of cell clusters were explored and cellular evolutionary processes were inferred by enrichment analysis and pseudotime analysis. In addition, a gene regulatory network (GRN) analysis was constructed to explore the regulatory role of transcription factors in AM progression. Results: We obtained the expression profiles of 42260 cells and identified 10 cell clusters. By comparing the differences in cell components between AM patients and controls, we found that significant abundance of endometrial cells (EC), epithelial cells (Ep), endothelial cells (En), and smooth muscle cells (SMC) in AM patients. Cell clusters with high CNV levels possessing tumour-like features existed in the ectopic endometrium samples. Moreover, the Ep clusters were significantly involved in leukocyte transendothelial cell migration and apoptosis, suggesting an association with cell apoptosis and migration. En clusters were mainly involved in pathways in cancer and apoptosis, indicating that En has certain malignant features. Conclusion: This study identified cell clusters with immune-related features, investigated the changes in the immune ecology of the microenvironment of these cells during AM, and provided a new strategy for the treatment of AM.

Background: Adenomyosis (AM) is a common benign uterine disease that threatens the normal life of patients. Cells associated with microenvironmental immune ecology are crucial in AM, although they are not as well understood at the cellular level.
Methods: Single-cell sequencing (scRNA-seq) data were used to construct an AM global single-cell map, to further identify relevant cell clusters and infer chromosomal copy number variation (CNV) in AM samples. The biological functions of cell clusters were explored and cellular evolutionary processes were inferred by enrichment analysis and pseudotime analysis. In addition, a gene regulatory network (GRN) analysis was constructed to explore the regulatory role of transcription factors in AM progression.
Results: We obtained the expression profiles of 42260 cells and identified 10 cell clusters. By comparing the differences in cell components between AM patients and controls, we found that significant abundance of endometrial cells (EC), epithelial cells (Ep), endothelial cells (En), and smooth muscle cells (SMC) in AM patients. Cell clusters with high CNV levels possessing tumour-like features existed in the ectopic endometrium samples. Moreover, the Ep clusters were significantly involved in leukocyte transendothelial cell migration and apoptosis, suggesting an association with cell apoptosis and migration. En clusters were mainly involved in pathways in cancer and apoptosis, indicating that En has certain malignant features.
Conclusion: This study identified cell clusters with immune-related features, investigated the changes in the immune ecology of the microenvironment of these cells during AM, and provided a new strategy for the treatment of AM. KEYWORDS adenomyosis, single-cell RNA sequencing, malignant cells, immune microenvironment, biological function Introduction Adenomyosis (AM) is a common benign gynecological syndrome, characterized by infiltration of endometrial glands and stroma in the myometrium (Bird et al., 1972). The prevalence of AM ranges from 5% to 70% (Graziano et al., 2015), with an average of 20%-35% of women worldwide suffering from AM (Abbott, 2017). The most common manifestations of AM are dysmenorrhea, infertility, and abnormal uterine bleeding (AUB), but some women with AM are asymptomatic (Peric and Fraser, 2006). From the epidemiologic data, AM can increase the risk of cancers, including ovarian, endometrial, breast, colorectal, and other cancers of women (Borgfeldt and Andolf, 2004;Kok et al., 2015). Furthermore, Bergeron previously reported that the definitive diagnosis of AM was based on the presence of ectopic endometrial tissue in the myometrium (Bergeron et al., 2006), but is now diagnosed by non-invasive techniques such as pelvic imaging (Chapron et al., 2020). Studies have also reported that the standard method of managing the disease is hysterectomy, while most patients desire to preserve their fertility (Stratopoulou et al., 2021). Despite improvements in diagnostic tools, awareness of AM remains poor (Leyendecker et al., 2015).
At the cellular level, the microenvironmental immune cells of AM play an important role. The number of macrophages, natural killer cells, and T cells in the endometrial stroma of AM increased significantly compared with women with mild focal AM or without the disease (Yang et al., 2004;Tremellen and Russell, 2012). Several malignant features exist in the epithelial cells (Ep) of AM, such as high migration capacity, which contribute to disease progression (Liu et al., 2021). Studies have confirmed that the endothelial cells (En) are damaged, and the uterus occurs the symptoms of bleeding, which is also important in adhesion and migration (Kruger-Genge et al., 2019). Furthermore, smooth muscle cells (SMC) have the ability to shrink and diastole, lack of contraction can cause uterine bleeding, which may lead to the occurrence of inflammation (Owens et al., 2004). It may thus be possible to comprehend the emergence of AM by concentrating on the mechanisms of change in cells related to the immunological microenvironment.
Single-cell RNA sequencing (scRNA-seq), an indispensable technique to dissect cellular heterogeneity and analyze cell types, can assist us in thoroughly comprehending the biological roles (Hedlund and Deng, 2018). Numerous effective methods to examine molecular alterations at the cellular level are provided by the scRNAseq (Tang et al., 2009). Moreover, research has demonstrated that rare clusters of AM were identified by scRNA-seq, confirming that the occurrence originates from endometrial migration (Liu et al., 2021). However, more studies are needed for further validation. In this study, we explored the states and transitions of the immune microenvironment cells of AM from a single-cell perspective. A comprehensive map of the AM single-cell ecosystem was depicted, relevant cell clusters were identified, and chromosomal copy number variation (CNV) was inferred for each AM sample. Furthermore, it was further confirmed the associated cluster of markers was involved in the signaling pathways and gene regulatory networks (GRN), which contributes to our understanding of the functions of the cluster markers in AM and at the cellular level.

Data sources
The AM scRNA-seq data including SRR12791871, SRR12791872, and SRR12791873 (Liu et al., 2021) were obtained from the Sequence Read Archive (SRA) of the National Center for Biotechnology Information (NCBI). A 50year-old woman with uterine fibroids, excluding the AM, and the endometrium tissue from this patient were used as a control sample. Moreover, two endometrium tissue samples were obtained from a 46-year-old AM patient and the samples were taken from eutopic endometrium (AM_EM) and ectopic endometrial (AM_EC) tissues.

Data preprocessing and construction of the single-cell atlas
We used the IntegrateData function (Butler et al., 2018) in the Seurat package (Stuart et al., 2019) to merge the scRNA-seq data, and performed cell clustering analysis according to default parameters. Uniform Manifold Approximation and Projection (UMAP) algorithm (Becht et al., 2018) was adopted for dimensionality reduction and visualization and mapped into single cell profiling. Subsequently, the FindAllMarkers function in Seurat package identified the specific marker genes in each cell cluster. Furthermore, the cell types underwent an immune response based on annotation and re-clustering of known marker genes.

Differential gene expression analysis
Differential expression analysis was performed based on the FindMarkers function in the Seurat package (Butler et al., 2018). Differentially expressed genes (DEGs) of different clusters in the Control, AM_EM, or AM_EC groups were identified. DEGs between normal tissues and AM tissues were screened by an adjusted p value <0.05 and |log fold change (log FC)|>0.5 being considered significant.

Evaluation of CNV in single cells
CNVs are primarily used to identify subclones of diseased cells and to infer tumor evolution (Yates and Campbell, 2012). The CNVs of each cell were assessed from the AM patients by the inferCNV package (inferCNV of the Trinity CTAT Project; https://github.com/broadinstitute/inferCNV) (Patel et al., 2014). To calculate the CNVs of AM_EM and AM_EC cells, the average or normal expression of genes from immune cells was applied as a reference and then determine the expression.

Functional enrichment analysis
To explore the biological functions involved in each cell cluster. We performed the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) on the clustered markers using the clusterProfiler package (Yu et al., 2012). p < 0.05 was considered statistically significant.

Pseudotime analysis
Pseudotime analysis could determine the dynamics of gene expression within cell types and trajectories over time (Cao et al., 2019), and infer cell evolution during the AM. A "branch" trajectory was constructed based on the Monocle three package (Trapnell et al., 2014) to explore the dysregulated changes in immune cells of AM patients and project the cells into low-dimensional space by UMAP, the parameters of Monocle three package are set to the default.

Construction of GRN
In this study, we constructed a GRN with transcription factors as the core to infer co-expression modules, to further explore the dysregulate regulatory mechanism of immune cellogenesis. Single cell regulatory network inference and clustering (SCENIC) (Aibar et al., 2017;Van de Sande et al., 2020) was used to infer gene regulatory networks based on singlecell expression profiles and identify cell states, providing important biological insights into the mechanism driving cell heterogeneity. Among them, the binding motifs of the transcription factors in the co-expression module were obtained from the JASPAR database (https://jaspar.genereg.net).

Statistical analysis
Statistical analyses were performed using R (https://www.rproject.org/). Expression levels of genes were analyzed using unpaired t-tests. If the p < 0.05 that considered statistical significance. The analyses in this study were based on the Bioinforcloud platform (http://www.bioinforcloud.org.cn).

Global single-cell atlas of the adenomyosis
To investigate the early cell population dynamics in AM patients, we analyzed the scRNA-seq data from endometrial tissue samples of AM patients and control donors. The analysis flow of this study is shown in Figure 1A, where we constructed a global single-cell landscape of AM. By cluster analysis, we divided 42,260 cells into 31 cell clusters, which were identified into 11 cell types based on known markers ( Figure 1B) (Supplementary Table S1), including endometrial cells (EC), fibroblasts, epithelial cells (Ep), endothelial cells (En), CD8 + T, CD4 + T, Naive T, macrophages (Mac), plasmacytoid predendritic cells (pDC), smooth muscle cells (SMC), and innate lymphoid cells (ILC). Among them, each cell marker exhibited a specific expression for the cell cluster ( Figure 1C). Chromosomal CNV analysis based on expression patterns at genomic intervals showed the presence of multicopy events in AM in ectopic endometrial samples ( Figure 1D). Further comparing the differences in cell composition between control and AM patients, we found that the highest abundance of En and EC was found in the AM_EC and AM_EM groups. However, the fibroblasts in the control group had the highest abundance ( Figure 1E). In summary, we delineated the single-cell profiles of AM patients to reveal the differences in microenvironmental cell components in AM patients.
Ecological landscape of adenomyosisassociated endometrial cells clusters AM occurs mainly in endometrial tissue, which has a higher cell abundance in EC of AM patients, therefore, subsequent studies will focus on this cellular cluster. We obtained 10 EC clusters by cluster analysis (Figure 2A). As shown in Figure 2B, almost all of these cell clusters were present in different groups of AM patients. Further exploration of the abundance of the cellular cluster in AM patients revealed a significant increase in the proportion of EC_TIMP3 cell clusters and a significant decrease in the proportion of EC_ZFAND2A clusters ( Figure 2C). Markers for the different clusters of the EC were mapped to the single-cell atlas, including ZFAND2A, KPT17, TIMP3, SPARCL1, PLAAT3, SPINT2, SCGB2A1, RGS5, CXCL2, and COL1A2 ( Figure 2D). Furthermore, EC clusters may be associated with cell motility, which was closely associated with focal adhesion and leukocyte transendothelial migration and apoptosis ( Figure 2E). Based on the pseudotime trajectory Frontiers in Genetics frontiersin.org    Frontiers in Genetics frontiersin.org 07 analysis, it was inferred that the EC_ZFAND2A cluster served as the developmental starting point, and then differentiated into other cell clusters ( Figure 2F). Furthermore, we constructed the GRN and found that the GRN with TFs as pivots was organized into five modules ( Figure 2G), such as RXRG, ZEB1, MSX2, DLX5, ELF5, ARNTL, to regulate the specific gene expression ( Figure 2H). Above all, we identified EC clusters of AM patients, defined marker genes for their specific expression, and elucidated the GRN of evolved EC clusters.

Ecological landscape of adenomyosisassociated en clusters
The abundance of En was significant in the AM_EC group, from which we inferred that En may play a facilitating role in the disease course. Therefore, the cluster analysis of En yielded 10 cell clusters ( Figure 3A). Further mapping of these cell clusters to the AM_EC, AM_EM, and control groups, and we found that mainly originated from AM patients ( Figure 3B). En_HLA-DRB1 and En_ID1 were increased in AM_EM, En_TPM1 was significant in the AM_EC group, and En_IFIT1 was significantly decreased in AM patients ( Figure 3C). Gene expression of markers for En clusters weas mapped to a single-cell atlas, including APOE, CLDN3, MMP1, HLA-DRB1, TPM1, ID1, IFIT1, ESM1, DES, and LUM ( Figure 3D). En clusters were mainly related to cell proliferation, such as pathways in cancer, apoptosis and extracellular matrix receptor interactions ( Figure 3E). Pseudotime trajectory analysis showed that the En_IFIT1 cluster evolved as a developmental starting point towards AM patients ( Figure 3F). The results of GRN for En clusters indicated that marker genes were divided into seven modules, and regulated by TFs, such as MYBL2, MAZ, and NEUROD2 ( Figure 3G). Figure 3H shows the transcription factors of En specific cell cluster. En cluster markers are regulated by transcription factors that promote the development of AM.

Exploring the ecological landscape in ep clusters of adenomyosis
Interestingly, Ep loses polarity and intercellular adhesion to gain migration capacity (Acloque et al., 2009), and the cell abundance of Ep clusters was significant in the AM_EC group. Therefore, cluster analysis of Ep clusters was again performed to obtain 10 cell clusters ( Figure 4A), which were mapped to AM_EC, AM_EM, and control groups according to their sample source ( Figure 4B). Compared with the control group, the cell abundance of Ep_ACTG2 was significantly increased in the AM_EM group, cell abundance of Ep_PALM2_AKAP2 was significant in the AM_EC group ( Figure 4C). Subsequently, we mapped the expression of cluster markers (PALM2_AKAP2, MEG3, ACTG2, LM07, WFDC2, MKNK2, S100A2, PMEL, CFD, and ESM1) to the single-cell atlas of Ep clusters ( Figure 4D). To further explore the biological signatures for the involvement of the Ep clusters in the AM, we performed enrichment analysis of the marker genes in the Ep clusters, showing that extracellular matrix receptor interactions, MAPK signaling pathway and apoptosis were significantly involved in Ep clusters ( Figure 4E). The developmental trajectory of Ep was explored by pseudotime trajectory analysis, and the results indicated Ep_PMEL cluster was in an early developmental stage and evolved into AM_MKNK2, AM_ESM1, AM_LM07, and AM_S100A2 ( Figure 4F). We further performed a GRN analysis of Ep clusters, showing that the marker genes of the Ep clusters were divided into five modules regulated by the transcription factors, such as KLF4, FOXP4, NFIA, and ERG ( Figure 4G). Furthermore, we explored the expression of these TFs in specific Ep clusters and found that KLF4 was the most highly expressed in the cluster ( Figure 4H), suggesting that high expression of KLF4 may be associated with the development of AM.

SMC-associated cell clusters of adenomyosis
During AM, ecological components of SMC clusters were significantly observed, so we will investigate their microenvironmental immune properties. The SMC clusters were re-clustered to obtain 10 cell clusters ( Figure 5A) and mapped to different samples ( Figure 5B). The SMC_TP53BP2 cluster had significant cellular abundance in the AM_EC and AM_EM groups compared to the control group, while SMC_IFI6 had significant cellular abundance in the control group ( Figure 5C). Next, we showed the expression of marker genes was significantly changed in SMC clusters for clusters ( Figure 5D). The SMC clusters were significantly involved in pathways, such as cytokine receptor interaction, vascular smooth muscle contraction and apoptosis ( Figure 5E). We also further explored the developmental trajectory of the SMC, clarifying the evolution trajectory from the SMC_TP53BP2 cluster to the SMC_IFI6, SMC_VCAN, and the SMC_CXCL8 cluster ( Figure 5F). In addition, the GRN analysis with TFs as the fulcrum yielded five modules ( Figure F5G) with gene expression of specific SMC regulated by REL, MAZ, and ETV7 ( Figure 5H). Taken together, these results suggest that certain specific clusters are closely associated with vascular smooth muscle contraction, can lead to smooth muscle ischemia in AM patients, and may promote the development of AM.
Frontiers in Genetics frontiersin.org

Discussion
To date, most studies of orthoendometrium in women with AM have focused on the expression of a single gene or a limited number of genes (Benagiano and Brosens, 2012;Xiaoyu et al., 2014). However, few studies have revealed the uniqueness of each cell in the AM process at the individual cellular level. In this study, we analyzed the scRNA-seq data from the endometrial tissue of one patient with AM and one fibroid patient as controls, and explored changes in the cellular state and immune microenvironment of AM. Ectopic endometrium samples present high CNV levels, and are considered a potential factor in the development of the AM. The presence of certain specific cell clusters was associated with the progression of AM. In conclusion, it is possible to explore the global status of AM patients at the cellular level, providing new insights into the indepth study of AM.
The ecological compositions of EC in the AM_EC and AM_EM groups were significant compared to the controls, indicating that EC plays an important role in the AM microenvironment. EC clusters were significantly involved in cell motility-related pathways such as focal adhesion and leukocyte transendothelial migration by enrichment analysis. AM patients suffer from uterine bleeding, pelvic pain, or infertility in women due to endometrial adhesion and destruction of the endometrium (Lacheta, 2019;Bulun et al., 2021). In the GRN, RXRG, ZEB1, MSX2, DLX5, and ELF5 regulated the EC clusters. Knockdown of SKP2 was found to reduce ZEB1 expression in endometrial stromal cells, thereby inhibiting their invasion and proliferation (Guo et al., 2022). However, few other TFs have been reported in AM. although many relevant mechanisms have been studied based on endometrial tissue, few studies have determined the effect of EC cells on AM at the cellular level. In conclusion, these results suggested that EC cells mediate and disrupt the endometrium, moreover, a cluster of markers regulated by TFs may promote the development of AM.
En has a dual role in immunology and pathology. On the one hand, the dysfunction will mediate the development of certain diseases, and on the other hand, they will actively mediate the immune response at the site of injury or infection (Pober and Sessa, 2007). Liu found co-localization of Ep and En markers in cluster one and promoted cell growth in AM (Liu et al., 2021), identical to the tumor-like characteristics reported by AM (Liu et al., 2018). Enrichment analysis showed that the En clusters were significantly involved in cancer-related pathways and extracellular matrix receptor interactions, indicating that the En was associated with cell proliferation and had certain malignant characteristics in the AM_EC group. Moreover, KLF4, FOXP4, NFIA, and ERG can regulate the markers of En clusters in the GRN. By inhibiting the biological functions of autophagy and metaphase during AM onset, KLF4 is abnormally reduced (Mei et al., 2022). However, other regulated TFs were almost rarely reported in AM, and whether En migrate needs further investigation.
Ep recognizes perturbations in their microenvironment, sends reinforcement signals, and transmits the signals to the immune system (Larsen et al., 2020). Studies have shown that epithelial immune cells of endometria can enhance cell survival and epithelial protective barrier function (Ho et al., 2006). Notably, the highest percentage of EP was found in the AM_EC and AM_EM groups; however, the number of EP was lower in the control group, which may be due to the significant postoperative endometrial thinning in the control group. Moreover, the Ep clusters were also significantly involved in the biological pathways such as leukocyte transendothelial migration, the MAPK pathway, and focal adhesions. Migration across the En after leukocyte adhesion, indicated that Ep was associated with high motility and migration. The MAPK pathway was required for the cell migration process, and the cytokines can also mediate cell migration (Liu et al., 2021). Furthermore, development of AM was improved by inhibiting the activated MAPK/ERK signaling pathway (Ying et al., 2021). In addition, study has shown that Ep loses their polarity and intercellular adhesion during adhesive epithelial interstitial transformation (EMT), to gain the ability to migrate to a mesenchymal phenotype (Acloque et al., 2009) and that EMT may play a key role in the pathogenesis of AM . In particular, during the conversion of Ep into En, a significant accumulation of angiogenic mimicry formation in AM_EC was found (Liu et al., 2021). Our results suggested that EP subsets significantly involved in pro-migratory pathways may play an important role in AM progression.
Abnormal proliferation of SMC in the endometriummyometrial junction area is an important cause of AM , and the emergence of AM causes hyperplasia and hypertrophy of the surrounding SMC (Zhai et al., 2020). Compared to the SMC in the normal uterus, the uterine SMC has hypertrophy and ultrastructural changes, which may have contractility effects on the myometrium (Mehasseb et al., 2010). In the AM global singlecell ecosystem, the components of the SMC are more prominent in the AM than in the controls due to AM-induced SMC abnormalities. SMC clusters were significantly involved in cytokine receptor interactions and vascular smooth muscle contraction. SMCs are capable of significant phenotypic changes in response to changes in local environmental cues, including cellcell and cell-matrix interactions, as well as various inflammatory mediators (Owens et al., 2004). In addition, vascular smooth muscle contraction occurs, leading to smooth muscle ischemia and dysmenorrhea in AM patients (Zhai et al., 2020). In GRN, we obtained five co-expression modules and three TFs, including REL, MAZ, and ETV7. Studies have demonstrated that REL was expressed and localized in the epithelial or stromal cells after castrated prostate patients (Rosa-Ribeiro et al., 2014). MAZ, as a transcriptional activator, may participate in the development of atherosclerosis (Ponnusamy et al., 2015). ETV7 promotes the resistance of breast cancer cells to chemotherapy and radiotherapy (Pezze et al., 2021). However, these transcription Frontiers in Genetics frontiersin.org factors have hardly been reported in AM. Generally speaking, maintaining the microenvironment homeostasis of SMC is critical for inhibition of development in AM.
In the present study, we found expression profiles of 42,260 cells and identified 10 cell clusters, including EC, fibroblasts, Ep, En, CD8 + T, CD4 + T, Naive T, Mac, pDC, SMC, and ILC. Among there, significant abundance of EC, Ep, En, and SMC in AM patients comparing the controls. Furthermore, the Ep clusters were mainly involved in leukocyte transendothelial cell migration and apoptosis; En clusters were mainly involved in pathways in cancer and apoptosis. Especially, some cell clusters were involved in cell migration and apoptosis may be promote the development of AM patients. Moreover, we obtained coexpression modules and TFs associated with the significant cell clusters by GRN comparing with the previous studies.
In conclusion, our study established a single-cell ecological landscape of the endometrium between control and AM patients, and explored the dynamic changes of immune cells during AM. However, this study has several limitations. First, the samples were too small and larger sample size is needed for a large-scale study. Moreover, this study was mainly based on bioinformatics analysis and therefore requires relevant molecular and cellular experimental validation.

Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/ Supplementary Material.