Immune system adaptation during gender-affirming testosterone treatment

Inclusion and ethics

Twenty-three adult individuals who were assigned female sex at birth and who were undergoing masculinizing gender-affirming treatment were enroled at specialist centres for transgender medicine in Stockholm, Uppsala, Linköping and Umeå in Sweden between 2016 and 2023. The study was approved by the Swedish Ethical Review Authority (2016/1422-31/1). Informed consent was obtained from all individuals. Only individuals who had not previously received testosterone treatment and who had normal sex hormone concentrations were included. Additionally, individuals with autoimmune diseases, immunodeficiencies or signs of continuing infection/inflammation were excluded from the study.

Testosterone therapy

Venous blood samples were collected at three timepoints: baseline, 3 months and 12 months following testosterone injections (Testosterone Undecaonate, Nebido administered once every 12 weeks). All patients received 1,000 mg of Nebido except for four participants who received reduced doses of 750 mg due to low body mass indices or haematocrit values. The baseline sample was collected within the 2 weeks before the start of testosterone treatment.

Measurement of serum sex hormones

Serum concentrations of sex hormones were analysed using liquid chromatography with tandem mass spectrometry assays at Gothenburg University as described previously⁶¹. The same method was used to analyse sex hormones from resulting culture supernatants as shown in Extended Data Fig. 4d. Briefly, calibrator stock solutions were prepared for all sex hormones and each internal standard stock solution was made separately using ¹³C₃-labelled versions of each steroid, except for dehydroepiandrosterone, which was labelled with d6. Steroid hormones were analysed using a two-dimensional liquid chromatography system consisting of an Acquity ultra-performance liquid chromatography system and a TQ-XS triple quadrupole mass spectrometer from Waters. The lower limit of detection (LLOD) was defined as the lowest peak with a signal more than three times the noise level. The lower limit of quantification (LLOQ) was defined as the lowest peak that was reproducible with a coefficient of variation of less than 20% and an accuracy of 80% to 120%. To circumvent problems with endogenous steroid amounts, the determinations of LLOD and LLOQ were performed in human serum pools with isotope-labelled steroids spiked at four different concentrations.

The amounts of bioavailable testosterone were calculated according to the formulae below⁶²:

in which rKb,ALB is the binding constant (0.601) for testosterone (TEST) to albumin (ALB), rKb,SHBG is the binding constant (1.0) for testosterone to sex hormone-binding globulin (SHBG) and P − ALB is a fixed value of 42.

Sample processing

A 4 ml sample of blood was drawn in EDTA-containing sterile vacutainer tubes from each participant in the sex reassignment therapy cohort and prepared as follows: 0.5 ml of blood was mixed with an equal amount of whole blood stabilizer⁶³ (Cytodelics AB), incubated for 10 min at ambient temperature and stored at −80 °C. A 1 ml aliquot of blood was mixed with PAXgene solution (BD Biosciences), incubated for 2 h at ambient temperature and stored at −80 °C. The remaining blood was centrifuged at 4 °C and 1,200g for 10 min, after which plasma was collected and stored at −80 °C. The leftover blood after plasma removal was mixed equally with PBS and layered over Lymphoprep (STEMCELL Technologies) for PBMC isolation by density gradient centrifugation following the manufacturer’s protocol. Cells were washed, counted and cryopreserved in a solution of 90% FBS (Sigma-Aldrich) mixed with 10% dimethylsulfoxide (DMSO; Sigma-Aldrich), initially stored at −80 °C overnight and then transferred to −150 °C for future use.

Bulk RNA-seq of whole blood samples

To analyse changes in gene expression, we performed RNA-seq using RNA extracted from PAXgene blood samples. The RNA samples were prepared using a QIAcube with the PAXgene Blood RNA Kit (Qiagen). Before cDNA library preparation, the quality of the RNA was assessed by determining the RNA integrity number using the Agilent 2100 Bioanalyzer with the RNA 6000 Pico Kit. The RNA concentration was measured using the Qubit Fluorometer with the Qubit dsDNA HS Kit (ThermoFisher Scientific).

For final sequencing and cDNA library preparation, an Advanta RNA-Seq XT NGS Library Preparation Kit was used with the Juno system (Standard BioTools Inc.). We performed Bulk RNA-seq on a NovaSeq 6000 instrument using one flow cell SP-200 (Illumina) with paired-end reads and a read length 2 × 100.

Data analysis of bulk mRNA-seq data

Bulk RNA-seq results from 59 samples from 20 individuals undergoing testosterone treatment were preprocessed with Kallisto⁶⁴. Quality control was provided by the National Genomics Infrastructure at Science for Life Laboratory, Stockholm, Sweden. To generate abundance estimates for all samples, the Kallisto program (v.0.46.2) was used to quantify abundances of transcript sequences in FASTA format using the Ensembl transcriptome Homo_sapiens.GRCh38.cdna.all.index (https://ftp.ensembl.org/pub/release-109/fasta/homo_sapiens/cdna/) for the Kallisto index. The Kallisto outputs were then imported into R using the tximport package, and the effect of ‘visit’ on whole blood mRNA expression was assessed using DESeq2 (ref. ⁶⁵) while accounting for interindividual variability and age effects. Before assessing differential gene expression, genes with fewer than 100 reads across samples were filtered out, as well as genes that did not have a normalized count of ten in at least one-fourth of the samples. The results from the differential gene expression analysis were used for gene set enrichment analysis of Hallmark pathways using clusterProfiler⁶⁶.

scRNA-seq experiments

Cryopreserved PBMCs obtained at baseline and after 3 months of testosterone treatment were thawed in thawing medium (RPMI 1640 HyClone supplemented with 10% FBS, 1% penicillin-streptomycin and Benzonase-nuclease (Sigma-Aldrich)). Cells were counted using a Cellaca MX (Nexcelom), plated and incubated for 1 h at 37 °C and 5% CO₂ to rest. Samples were then either left untreated or stimulated ex vivo with LPS (100 ng ml^−1) or R848 (1 μg ml^−1) for 4 h. After stimulation, the cells were collected, and supernatants were stored for later analysis by SIMOA (Quanterix)⁶⁷.

Viability and cell counts were assessed after resuspending collected cells in PBS with 0.04% BSA (ThermoFisher Scientific). The cells were then prepared for scRNA-seq using the 10x Genomics 3′ v.3.1 (dual index) kit according to the manufacturer’s instructions (catalogue no. CG000315 Rev B) on a Chromium Controller. Approximately 1 × 10⁴ cells from each condition were loaded onto separate wells of a 10x Genomics chip and the Chromium Controller was used to create GEM emulsions. The target recovery was 6,000–7,000 cells per condition. The libraries were sequenced on an Illumina NovaSeq 6000 platform, using paired- end reads (configuration 28 × 10 × 10 × 90) with 20,000 reads per cell.

scRNA-seq data analyses

CellRanger with default parameters was used to process FASTQ-files and align sequencing reads from 10x Genomics 3′ HT v.3. and 3′ GE towards the human genome. Cells were further filtered using a bimodal distribution-based approach, excluding those with read counts below (considered low quality) or above (considered technical artifacts) cut-off thresholds. The cut-off thresholds for each sample were chosen on the basis of distribution shape of read counts to retain biologically relevant cells and to eliminate technical artifacts. Cells with mitochondrial gene expression above 10% were also filtered out. All scRNA-seq data were preprocessed in Python using Scanpy v.1.9.1. For each sample, normalization by counts per cell (target sum = 1 × 10⁴) and feature scaling were applied to the CellRanger outputs for each sample, followed by linear dimensionality reduction using PCA and uniform manifold approximation and projection (on top 2000 variable genes), nearest neighbours (n = 10) computation and identification of clusters (res = 1). Clusters were annotated on the basis of canonical marker genes. BTMs⁶⁸ were used to compare transcriptional patterns before and during testosterone treatment and in response to stimulation.

NicheNet analyses

The NicheNet analysis and circus plots were created following the standard workflow available from NicheNet⁶⁹ and circlize⁷⁰. Specifically, differentially expressed genes between samples from baseline and after 3 months of testosterone treatment were identified using Seurat’s (v.4.3.0) built in function FindMarkers and filtered with an adjusted P value of less than 0.05 and an absolute value for the average fold change of at least 0.15. Ligand activities were calculated, and the top upstream ligands that could explain the observed target gene expression changes were selected. The ligand–target links were filtered on the basis of their weights (strength of the ligand–target relationship), with links belonging to the lowest 66% of scores being removed. The circos plot blocks were coloured according to a gene’s target cell, inferred as the cell type with the highest mean-value change between the two visits. The widths of the blocks indicate the potential of each receptor to be influenced by all shown ligands, with some interactions not visible due to the cut-off weight threshold. The transparency of the arrows indicates the regulatory potential of a ligand–target interaction (the more transparent, the weaker the regulatory potential).

sc-ATAC-seq and data analysis

sc-ATAC-seq experiments were conducted on the 10x Chromium platform, following a previously described protocol⁷¹. Briefly, cells were washed with PBS containing 0.04% BSA and nuclei subjected to isolation as per the manufacturer’s instructions. After counting, approximately 10,000 nuclei were used for tagmentation. The tagmented nuclei were then loaded for capture using the 10x Chromium controller. Following gel emulsion generation, we carried out linear amplification and DNA purification according to the manufacturer’s protocol. The resulting DNA was used for library construction, following the guidelines provided on the manufacturer’s website. The libraries were quantified using an Agilent Bioanalyzer and sequenced on an Illumina NovaSeq S4 sequencer, with the following setup: 50 bp read 1N, 8 bp i7 index, 16 bp i5 index and 50 bp read 2N. In this setup, 1N and 2N refer to the DNA insert sequencing, while i5 and i7 sequencing identify the individual barcodes of single cells.

The 10X Genomics cellranger pipeline (cellranger-atac mkfastq, count and aggr) was followed for 10x sc-ATAC-seq analysis. Cellranger aggr outputs were used for downstream analysis in R using the Signac package. We performed quality control using Signac’s default settings for transcriptional start site enrichment score, nucleosome banding pattern, sequencing depth and complexity, and fraction of fragments in peaks. The ratio of reads in genomic blacklist regions was calculated using the FractionCountsInRegion function with the blacklist for hg38. After quality control, a total of 143,624 peaks (features) across 12,773 cells remained for further analysis. The number of cells per sample varied between 636 and 4,632 for the eight total samples analysed. We applied frequency-inverse document frequency normalization, followed by feature selection and dimensionality reduction using singular value decomposition on the frequency-inverse document frequency matrix. We performed uniform manifold approximation and projection dimensionality reduction⁷² on the first 30 latent semantic indexing components, with latent semantic indexing components capturing technical variation excluded from further analysis. K-nearest neighbour graph construction and clustering were done using the smart local moving algorithm, resulting in the identification of 21 unique clusters. Gene activities were used for cluster annotation, with gene activities determined using the GeneActivity function followed by log normalization. Five main immune clusters were identified and used for further analyses. TF motif analysis was conducted by identifying overrepresented motifs in a set of differentially accessible peaks between pre- and post-testosterone therapy (3 or 12 months) for all the five immune subsets using hypergeometric tests and P values corrected for several hypotheses (Benjamini–Hochberg).

Immune cell profiling by mass cytometry

Blood samples were mixed with a stabilizer⁶³ (Whole blood processing kit component; Cytodelics AB) within the first hour post blood-draw and cryopreserved according to the manufacturer’s recommendations. Samples were then thawed, fixed and lysed using Lysis and Wash buffers (Whole blood processing kit; Cytodelics AB). After fixation/lysis, 1–2 × 10⁶ cells per sample were plated and cryopreserved using CRYO#20 (Cytodelics). For staining, cells were thawed at 37 °C, barcoded using an automated liquid handling robotic system (Agilent Technologies)⁷³ using the Cell-ID 20-plex Barcoding kit (Standard BioTools Inc.) as per the manufacturer’s recommendations and stained batch-wise after pooling. Cells were washed using cell staining buffer (CSB) (Standard BioTools Inc.), FcR blocked using an in-house-prepared blocking solution for 12 min at ambient temperature then stained using a cocktail of metal-conjugated antibodies targeting surface antigens (Broad extended panel) and incubated for 30 min at 4 °C. Cells were washed twice with CSB and fixed overnight using 2% formaldehyde in PBS (VWR international). The panel of antibodies used is listed in Supplementary Table 1.

For cells from whole blood pretreated and stimulated in vitro, we performed intracellular staining. Cells were first stained with a cocktail of antibodies targeting surface antigens (Supplementary Table 2) and then washed twice with CSB, fixed and permeabilized using Foxp3/Transcription Factor Staining Buffer Set (ThermoFisher Scientific) according to the manufacturer’s instructions. Cells were then stained using a cocktail of metal-conjugated antibodies targeting intracellular antigens (Supplementary Table 3) and incubated for 1 h at ambient temperature. Cells were washed twice with CSB and fixed overnight using 2% formaldehyde in PBS.

For acquisition by CyTOF XT⁷³, cells were stained with DNA intercalator (0.125 mM Iridium-191/-193 or MaxPar Intercalator-Ir (Standard BioTools Inc.) in 2% formaldehyde and incubated for 20 min at ambient temperature. Cells were washed twice with CSB followed by two washes with Maxpar Cell Acquisition Solution (CAS) Plus (Standard BioTools Inc.) before being filtered through a 35 mm nylon mesh, diluted to 500,000 cells ml^−1 using CAS Plus and divided into polypropylene tubes. A total of 2 × 10⁶ cells per tube in pelleted form were then placed in the chilled carousel of the CyTOF XT instrument (Standard BioTools Inc.). EQ Six (EQ6) element calibration beads (Standard BioTools Inc.) were added to a tube and placed in the carousel. The autosampler of the CyTOF XT dispensed CAS Plus to the pelleted sample tubes, mixed with EQ beads 0.1×, and then acquired on CyTOF XT mass cytometers at a rate of 300–500 cells s^−1 using CyTOF software v.8.0 with noise reduction, event length limits of 10–150 pushes, and a flow rate of 0.030 ml min^−1.

Mass cytometry antibodies and reagents

Purified antibodies were obtained in carrier/protein-free buffer and coupled to lanthanide metals using the MaxPar X8 or MCP9 antibody conjugation kits (Standard BioTools Inc.) as per the manufacturer’s recommendations. Metal-conjugated antibodies were also purchased from Standard BioTools. The antibodies used for this study are listed in Supplementary Tables 1–3.

Mass cytometry data analyses

Samples from participants undergoing sex reassignment therapy were processed through mass cytometry in four batches to investigate immune composition and phenotype. This involved analyzing.fcs files from 60 samples from 20 series of participants receiving testosterone treatment. Data analysis was conducted in R. The data were arcsin h transformed with a cofactor of five using the flowCore package. Beads and dead cells were filtered out. Batches were combined, and batch effects in marker expression were eliminated using the sva package. The resulting matrix was used for immune composition analysis with the FlowSOM package⁷⁴.

Initially, 30 clusters were identified, neutrophil clusters were annotated, and the remaining non-neutrophil cells were clustered into a total of 100 clusters. A total of 113 unique clusters were annotated on the basis of median marker expression using the pheatmap package. A total of 12,377,068 cells from the 60 samples of participants undergoing testosterone treatment were further analysed. This analysis included investigating immune phenotypes using PAGA⁷⁵ (see below) and examining the effects of testosterone on immune cell composition using a mixed-effects model with the lme4 package. For linear mixed-effects models, the frequency of 35 immune subsets was modelled considering visit (baseline, 3 months and 12 months) and age as fixed effects, and participant ID as random effect. Significant visit effects were determined using a P value of 0.05 and a 5% FDR threshold, with beta coefficients indicating the directionality of the effect.

Spectral flow cytometry analysis of AR and ESR expression

For ESRa staining, PBMCs were extracted from heparinized whole blood, as described above. One million live cells were aliquoted per test, washed twice in ice-cold PBS, and incubated with LIVE/DEAD Fixable Blue dye (ThermoFisher Scientific) for 10 min at 4 °C. PBMCs were then washed in ice-cold FACS buffer (2% FBS, 0.5 mM EDTA in PBS) and FcR blocked using an in-house-prepared solution for 10 min at ambient temperature. The Horizon Brilliant Stain Buffer Plus (BD Biosciences) and extracellular antibodies (Supplementary Table 4) were added, and cells were incubated for 30 min at 4 °C followed by fixation and permeabilization with Fixative buffer (Cytodelics AB) or Foxp3/Transcription Factor Staining Buffer Set (ThermoFisher Scientific) according to the manufacturer’s instructions. Cells were mixed with FcR block buffer and, after 10 min at room temperature, intracellular antibodies (Supplementary Table 4) were added, and the samples were incubated for 30 min at 4 °C. For AR staining, WBCs from heparinized whole blood were prepared using a Cytodelics kit, and 1.5 million fixed-permeabilized cells was aliquoted per test and exposed to FcR block (BD Biosciences) for 30 min at 4 °C. Horizon Brilliant Stain Buffer Plus (BD Biosciences) and all-antibody cocktail (Supplementary Table 4) were added, and cells were incubated overnight at 4 °C. AR and ESRa antibody concentrations were established on the cell line MCF7 (ATCC); the specificity of AR antibody was also verified using a competitive displacement approach on MCF7 cells. Briefly, 60,000 cells were collected at passage two, fixed-permeabilized using Cytodelics kit, FcR blocked and stained as described for WBCs. Unconjugated antibodies and isotype controls information is present in Supplementary Table 4. After a wash in cold FACS buffer, data were acquired using an Aurora spectral cytometer (Cytek Biosciences). Cytobank Community (Beckman Coulter) software was used for data analysis.

PBMC stimulation and intracellular staining by spectral flow cytometry

Cryopreserved PBMCs obtained from individuals undergoing gender-affirming testosterone treatment were collected at baseline and after 3 months of testosterone treatment. These cells were thawed in thawing medium (RPMI 1640 HyClone supplemented with 10% FBS, 1% penicillin-streptomycin and Benzonase-nuclease (Sigma-Aldrich).

The cells were then counted using a Cellaca MX (Nexcelom), plated and incubated for 1 h at 37 °C and 5% CO₂ to rest. After this, some samples were left untreated while others were stimulated ex vivo with PMA (50 ng ml^−1) and Ionomycin (1 μg ml^−1) for 4 h. Brefeldin A (5 μg ml^−1) and Monensin (2 μg ml^−1) were added during the last 3 h of stimulation.

Following stimulation, the cells were washed twice in ice-cold PBS and then incubated with LIVE/DEAD Fixable Blue dye for 10 min at 4 °C. The cells were then washed in ice-cold FACS buffer and FcR blocked using blocking solution prepared in-house for 10 min at ambient temperature.

Horizon Brilliant Stain Buffer Plus was added, and the cells were stained with a cocktail of fluorochrome conjugated antibodies targeting surface antigens for 30 min at 4 °C (Supplementary Table 4). The cells were then fixed using Fix, Lysis and Wash buffers (Whole blood processing kit; Cytodelics AB) and permeabilized using permeabilization buffer (ThermoFisher Scientific).

Next, the cells were stained with a cocktail of antibodies targeting intracellular antigens (Supplementary Table 4) for 30 min at 4 °C and then acquired using an Aurora spectral cytometer.

Plasma protein profiling by Olink

Plasma protein data was generated using the Olink assay, a proximity extension assay (Olink AB)⁷⁶. Plasma (20 μl) from each sample was thawed and analysed using a Target Inflammation panel (Olink AB), at the Affinity Proteomics Stockholm, Science for Life Laboratory or Olink AB. In these assays, plasma proteins are dually recognized by pairs of antibodies coupled to a cDNA-strand that ligates when brought into proximity by its target, extended by a polymerase and detected using a Biomark HD 96.96 dynamic PCR array (Standard BioTools Inc.). Analyses of differentially abundant plasma proteins were performed using linear mixed-effects models considering age as fixed effects.

Whole blood pretreatment in vitro using testosterone and AR antagonist for Olink analysis

A blood sample obtained from a healthy female donor was mixed in equal ratio with WB-STIM buffer (Cytodelics AB) without phenol red. The sample was then split into three groups: untreated, treated with testosterone (Sigma-Aldrich) alone at 10 ng ml^−1, or treated with a combination of testosterone and the AR antagonist Enzalutamide (Sigma-Aldrich) at 2.3 μg ml^−1. All samples were incubated for 28 h at 37 °C and 5% CO₂. After incubation, supernatants were collected, cryopreserved and later analysed using the Olink Target Inflammation panel (Olink AB) as described above.

Whole blood pretreatment and stimulation in vitro for Nanostring and mass cytometry analysis

For the in vitro pretreatment step, blood samples were mixed in equal ratio with WB-STIM buffer (Cytodelics AB) without phenol red and split as follows: untreated, treated with DHT (Sigma-Aldrich) alone at 10 ng ml^−1, treated with DHT combined with Enzalutamide (Sigma-Aldrich) at 2.3 μg ml^−1 or treated with fulvestrant (Sigma-Aldrich) alone at 100 nM. Samples were incubated for 20 h at 37 °C and 5% CO₂. DHT was chosen because this androgen cannot be converted to oestradiol by aromatase⁷⁷. Fulvestrant is a degrader of the ESR and blocks oestradiol-mediated signalling⁷⁸.

For Nanostring analyses, blood samples from healthy cisgender female donors (n = 11) were pretreated and then immediately stimulated with either LPS (10 ng ml^−1) or R848 (1 μM) for 3 h or left unstimulated as a control. Samples were then centrifuged at 4 °C and 1,200g for 10 min and supernatants were collected, cryopreserved and analysed using SIMOA. The remaining 1 ml of blood was mixed with PAXgene solution (BD Biosciences), incubated for 2 h at ambient temperature and stored at −80 °C. RNA samples were prepared using a QIAcube with the PAXgene blood RNA kit (Qiagen) and analysed using the Nanostring nCounter Sprint Profiler system with a broad human immune response panel (Human Immunology v.2 Gene Expression CodeSet) as described previously⁶. For each sample, 100 ng of total RNA in a final volume of 5 μl was mixed with a capture probe and a reporter probe tagged with a fluorescent barcode from the gene expression code set. Probes and target transcripts were hybridized overnight at 65 °C for around 19 h according to the manufacturer’s recommendations. Hybridized samples were run on the Nanostring nCounter instrument using the corresponding protocol, in which excess capture and reporter probes were removed and transcript-specific ternary complexes were immobilized on the surface of the cartridge. The images from samples were scanned at high resolution by the nCounter instrument and gene expression data were collected after scanning and image processing.

For mass cytometry analyses of cytokine production, blood samples from healthy cisgender females (n = 5) of reproductive age were collected before the ovulation phase of the menstrual cycle (day 1–10 from the first day of menstruation), pretreated and then immediately stimulated with either LPS (0.1 ng ml^−1) or PMA (50 ng ml^−1) combined with ionomycin (1 μg ml^−1) for 4 h or left unstimulated as a control. Brefeldin A (5 μg ml^−1) and Monensin (2 μg ml^−1) were added in all conditions. Samples were then fixed and lysed using Lysis and Wash buffers (Whole blood processing kit; Cytodelics AB). After fixation/lysis, cells were cryopreserved using CRYO#20 (Cytodelics AB) and analysed using intracellular staining mass cytometry as described above.

Analyses of Nanostring gene expression data

Batch-normalized data were log-transformed and scaled to have unit variance and zero mean. This was followed by principal component analysis (PCA). The resulting PCAs were then plotted alongside the PCA loadings of hallmark TNF genes.

Quantification of IFNa and IFNb by Simoa

IFNa subtypes were quantified in plasma and in supernatants of ex vivo–stimulated PBMCs using Simoa digital ELISA (Quanterix) with HomeBrew assays as previously described⁷⁹. Several IFNα subtypes were measured using a pan-IFNα subtype assay (Quanterix), with IFNa17 (PBL Assay Science) as a reference standard. Antibodies cloned from two patients with mutated APS1 were used for multi-IFNα subtype quantification. The 8H1 clone was coated on paramagnetic beads and used as the capture antibody (0.1 μg ml^−1), and the 12H5 clone was biotinylated at a ratio of 30:1 and used as the detector. The limit of detection for IFNα was 0.03 fg ml^−1. IFNβ was also quantified in plasma from the cohort. For the IFNβ assay, the 710906-9 IFNβ antibody (PBL Assay Science) was coated on paramagnetic beads (0.3 μg ml^−1) and used as a capture antibody. The 710323-9 antibody (PBL Assay Science) was biotinylated and used as the detector (30:1). Recombinant IFNβ (PBL Assay Science) served as a standard to determine unknown sample concentrations. The LOD for IFNβ was 0.3 pg ml^−1.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.