The data used in this paper were obtained from two main sources, the TCGA and GTEx databases. All transcript RNA-seq datasets of tissues were downloaded from the University of California Santa Cruz (UCSC) Xena database (30). Transcript expression levels were estimated using the FPKM (fragments per kilobase per million fragments mapped) method in HTSeq. A total of 320 samples (130 tumor tissues and 190 normal tissues) were included in our study. All specimens were taken from the thyroid. Data on 130 tumor tissues and 17 normal tissues were obtained from the TCGA database, and data on 173 normal tissues were obtained from the GTEx database. Clinical and follow-up data for each tumor sample (n = 130) were obtained from the TCGA data portal (up to Jan 28, 2016). The data are summarized in Supplementary Table 1 .

Batch correction was performed on the RNA-seq datasets from different sources. Batch correction was performed using the Combat method of the sva package (version 3.36.0) in R. Before batch correction, the normalization of the GTEx data was adjusted from log2(x+0.001) to log2(x+1) to align with the normalization of the TCGA data.

The limma package (version 3.44.1) was used to identify DEGs (31). |Fold Change (FC)| > 2 and p-value < 0.05 after FDR (false discovery rate) adjustment were established as cutoff criteria. Heatmaps of DEGs were generated using the R package pheatmap (version 1.0.12).

Weighted gene coexpression network analysis (WCGNA) was used to construct a coexpression network based on scale-free topology to detect coexpression modules of genes with similar patterns of expression (32). In our study, DEGs were input into the WGCNA package (v1.69) in R to construct weighted coexpression modules. The module membership (kME) value was used to weight the importance of a gene in a module (33). In our study, DEGs with |kME| values > 0.8 were considered hub genes of the weighted gene coexpression network (34). The relationships between modules and clinical traits were confirmed by Pearson’s correlation and are shown in a heatmap. The following clinical traits were analyzed: T, tumor stage, divided into two groups according to the 8th American Joint Committee on Cancer (AJCC) editions (T1–2 and T3–4); N, lymph node metastasis stage, divided into two groups according to the 8th AJCC editions (N0/Nx and N1); M, distant metastasis stage, divided into two groups according to the 8th AJCC editions: (M0/Mx and M1); stage, cancer stage, divided into two groups according to the 8th AJCC editions (stage I–II and stage III–IV); age, divided into two groups (≤60 years of age and >60 years of age); and extent of the primary tumor. In thyroid surgical practice, the thyroid is divided into three anatomical zones, the left lobe, right lobe, and thyroid isthmus. Patients were divided into two groups according to postoperative pathological findings depending on whether they had tumors confined to a single lobe or tumors in multiple lobes. A p-value < 0.05 was considered significant.

Protein–protein interaction (PPI) analysis is useful for tracking potential interactions among DEGs. We used the STRING online database (version 11.0) to construct a PPI network of the DEGs (35, 36). Degree indicates the number of connections/interactions involving the hub genes. To screen the hub genes of the PPI network, the selection criterion for nodes was set as degree ≥10 (37).

The intersection of hub genes of the weighted gene coexpression network and PPI network were regarded as hub genes of TC in males. A network of the 47 overlapping hub genes and 13 modules containing these hub genes was constructed by Cytoscape (version 3.8.0). Genes within the same module were analyzed for enrichment of Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways using DAVID bioinformatics resources (version 6.8) (38, 39). p-values <0.05 and counts >3 were considered threshold values (40).

Survival analyses were performed using the survival package in R (survival package version 3.1-8). Patients were split into two groups (low expression/high expression) according to the expression level of each hub gene. The Kaplan–Meier method and a log-rank test were used to compare the DFS rates between the two groups. Additionally, univariate Cox proportional hazard regression analysis was performed for each hub gene to test its association with the DFS rate of male TC patients. Hub genes with p-values less than 0.05 in the Kaplan–Meier survival curve and univariate Cox analysis were entered into a multivariate Cox proportional hazards model. A p-value <0.05 was set as the threshold value. Finally, ASF1B expression was screened using the method described above. Then, a ROC curve was used to establish the best cutoff point for ASF1B levels that optimally predicted the DFS rate. The best cutoff value for the predictive equation was achieved when the Youden index (YI = sensitivity + specificity − 1) was at its maximum. ASF1B high/low expression groups were determined according to the best cutoff value. Finally, the predictive power of ASF1B expression levels and clinical features in combination was analyzed using multivariate Cox regression. p < 0.05 was used as the criterion for statistical significance. The area under the curve (AUC) values and ROC curves were calculated using SIGMAPLOT (version 14, Systat Software INC, CA, USA).

To analyze the relationship between the tumor immune microenvironment and the expression levels of ASF1B, TISIDB online tools, the CIBERSORT algorithm, and the ESTIMATE algorithm were used. TISIDB is a database for tumor and immune system interactions (41). We evaluated the difference in the ASF1B expression level between responders and nonresponders to anti-PD-L1 (atezolizumab) therapy using TISIDB online tools. The data used for this analysis were obtained from a previously published urothelial cancer research study (42).

We further explored the relationship between the ASF1B expression level and tumor immune cell infiltration in TC in males using the CIBERSORT tool. CIBERSORT is a computational inference tool that can estimate the abundance of 22 leucocyte subtypes in tumors based on RNA-seq data. Each tumor sample was scored by the “CIBERSORT” R package (version 1.0.3).

In this study, correlations between ASF1B expression levels and infiltration level/ESTIMATE scores were calculated using the Spearman method. For consistency with the DEG results, the ASF1B expression levels were normalized to FPKM log2(x+1) format for the Spearman correlation analysis. The R package “ESTIMATE” (version 1.0.13) was used for this analysis, and a p-value <0.05 was set as the threshold value.

To verify the difference in the ASF1B expression level between male TC tissues and normal male thyroid tissues, we downloaded relevant data from the GEO database. The validation set consisted of 96 samples from male TC patients (71 tumor tissues and 25 normal tissues). Gene expression profiling data were taken from the GSE-29265, GSE-3467, and GSE-76039 datasets. These data were quantified using the metric log2(FPKM) and then corrected for batch effects using the Combat method of the R package sva (version 3.36.0). The differences in ASF1B expression levels between tumor tissues and normal tissues were compared by independent-samples t-tests in R software (version 3.6.3). A p-value <0.05 was set as the threshold. Violin plots were utilized to visualize the distribution of ASF1B expression levels among samples.

The human thyroid cell line (KTC-1) was purchased from Zolgene Technology Co., Ltd. (Zolgene, China). The KTC-1 cell line was originally derived from a 68-year-old white male TC patient (43). This cell line is consistent with the main characteristics of the samples in our study. KTC-1 cells were cultured in RPMI 1640 (Meilune Biotech, China) supplemented with 10% fetal bovine serum (FBS), 100 U/ml penicillin, and 100 μg/ml streptomycin. All cells were maintained in a humidified incubator with 5% CO₂ at 37°C.

Short hairpin RNA (shRNA) targeting ASF1B and scrambled control shRNA were designed and synthesized by Anburui Co., Ltd. (Anburui, China). All shRNA sequences used in this study are provided in Supplementary Table 6 . KTC-1 cells were seeded at 1 × 10⁴ cells/well into 6-well plates in complete media. The shRNAs and the shRNA-negative control vector were purchased from Shanghai Yusheng Biotechnology (Shanghai, China). These vectors were introduced into cells using Hieff Trans™ liposomal transfection reagent following the manufacturer’s protocol. The cells were washed and collected for subsequent experiments after incubation at 37°C for 12 h.

ASF1B knockdown levels were measured in KTC-1 cells transfected with 3 different shRNAs targeting ASF1B (sh1, sh2, and sh3). A blank group (Blank) and a negative group (shRNA-NC) were included as controls.

RNA was extracted from cells using RNAiso Plus (Takara, China). A total of 1 µg of RNA was transcribed to cDNA using NovoScript^® Plus All-in-one 1st Strand cDNA Synthesis SuperMix (gDNA Purge) Novoprotein, E047-01A (Shanghai, China). The resulting RT product was amplified using NovoStart^® SYBR qPCR Supermix (Novoprotein, China), and the data were collected using an iQ5 Real-Time PCR Detection System (Bio-Rad, USA). Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) was used as the loading control. The relative gene expression was generated using the 2^−ΔΔCT method (ΔCT target gene − ΔCT control gene).

Cellular proliferation was determined by CCK-8 following the manufacturer’s instructions (Meilune, China). Cell suspensions (100 µl, 5,000 cells) of sh-ASF1B and sh-ASF1B NC cells were added to each well of a 96-well plate. After cell culture for 0 h, 24 h, 48 h, and 72 h, respectively, 10 μl of CCK-8 solution was added to each well and incubated at 37°C for 1 h. The absorbance was measured at a wavelength of 450 nm for each well.

We analyzed the apoptosis and cell cycle of KTC-1 cell lines using Flow Cytometry. After indicated treatments, for the cell apoptosis analysis, KTC-1 cells were collected separately, and incubated with Annexin-V (Meilune, China) in the dark for 15 min, and the apoptotic cells were analyzed using a FACSCalibur (BD, USA). For the cell cycle analysis, KTC-1 cells were collected separately, resuspended in 1 ml of ice-cold ethanol 75% overnight at 4°C, incubated with RNaseA (Meilune, China) at 37°C for 30 min, incubated with propidium iodide (PI) in the dark at 4°C for 30 min, and then detected by FACSCalibur (BD, USA).

The Transwell invasion assay was performed using a Transwell chamber with a Matrigel-coated filter (Corning, USA). First, cells were suspended in serum-free medium. Then, the cells were added to the upper chamber and cultured with serum-free DMEM, whereas the lower chamber was filled with complete medium supplemented with 20% fetal bovine serum (Gibco, USA) and 1% penicillin and streptomycin (HyClone, UK). After being cultured in a 24-well plate at 37°C and 5% CO₂ for 24 h, the cells were blocked with 4% paraformaldehyde fixative and dyed with crystal violet (20% methanol, 0.1% crystal violet in H₂O). Finally, the invasive cells were counted in three random fields under a microscope. The experiment was repeated three times.

The KTC-1 cells were selected, digested and counted, and then inoculated into 6-well culture plates. Scratch was made by a 200-μl pipette tip when the cell confluence grew to 90%. Select the same point of view to take a picture at 0 h and 12 h. The wound width was monitored by an inverted microscope for assessing the capacity of KTC-1 cell migration.

Protein cells were collected using RIPA lysis buffer (KGP702, KeyGen Biotech, China). Protein concentrations were determined using a Pierce BCA Protein Assay Kit (Thermo, No. 23227, USA). SDS-PAGE (10%) was applied to isolate proteins, which were later transferred onto PVDF membranes. Membranes were blocked and incubated with primary antibodies overnight in nonfat milk at 4°C. Secondary antibody was applied at RT for 1–1.5 h. Signal detection was performed with chemiluminescence imaging system (Bio-Rad, USA).

A rescue experiment was performed by overexpressing FOXP3 in ASF1B-silenced KTC-1 cells. Cells were transfected with FOXP3-overexpressing plasmid (zolgene Biotech, China; Plasmid sequences are listed in Supplementary Table 7 ) or empty plasmid using Lipofectamine3000 kit (Invitrogen, USA) according to the manufacturer’s instruction. Cells were harvested at 48–96 h after transfection and RT-qPCR was used to detect mRNA expression. The effects of FOXP3 overexpression on proliferation and migration of ASF1B-silenced KTC-1 cell lines were examined by CCK8 assay and scratch assay; the detailed experimental scheme is the same as before.

Figure 1 shows the workflow of our research. Since the microarray data were collected from two different sources, data were first subjected to batch correction before differential gene expression analyses ( Supplementary Figure 1 ). A total of 2,933 genes were differentially expressed, with 331 upregulated genes and 2,602 downregulated genes ( Figure 2 and Supplementary Table 2 ). The heatmap represented the top 10 upregulated DEGs and top 10 downregulated DEGs between tumor tissues and normal tissues ( Figure 3 ).

A total of 130 samples with clinical data available and 2,933 DEGs were included in the WCGNA. Significant abnormal values were not found in the hierarchical clustering results ( Supplementary Figure 2A ). To ensure a scale-free network, soft-threshold analysis was performed on all samples. We chose a soft thresholding power (β) of 4, and this value fitted the scale-free topology assumptions. The linear regression model fitting index R2 was greater than 0.8. The mean connectivity was almost 0 when β = 4 ( Supplementary Figures 2B, C ). Using the dynamic tree cut and merging similar modules method, a total of 23 modules were identified ( Figure 4 ). A heatmap shows the relationships between the clinical traits and each module by Pearson’s correlation. The turquoise module (r = 0.45, p = 1e-6) and blue module (r = −0.44, p = 4e-6) were strongly correlated with lymphatic metastasis; blue modules (r = −0.27, p = 0.04) were strongly correlated with T stage ( Figure 5 ). |kME| value > 0.8 was used as a screening standard for the selection of hub genes and identified 222 hub genes among the DEGs ( Supplementary Table 3 ).

We constructed a PPI network of 2,933 DEGs using the STRING database. A total of 300 genes with a degree ≥10 were selected from the PPI network as hub genes ( Supplementary Table 4 ).

A total of 47 overlapping hub genes between the WGCNA and PPI networks were identified using the Venn diagrams online tool (44) ( Supplementary Figure 3 ). Forty-seven hub genes were distributed across 13 modules ( Figure 6 ). Gene ontology (GO) and KEGG enrichment analyses were conducted for the genes in each module, and the gene sets of 9 modules showed enrichment in GO or KEGG pathways ( Supplementary Figure 4 ). The most significant pathway for each module is indicated in the text label near each module in Figure 6 . Different colors represent different modules and the genes they contain.

Kaplan–Meier survival and univariate Cox analysis showed that a high ASF1B expression level was a risk factor for DFS (the patients were grouped by median ASF1B expression) ( Figure 7A and Supplementary Table 5 ). The best cutoff point of ASF1B expression (2.96) was calculated according to the maximum Youden index in the ROC curve for the DFS rate. The sensitivity and specificity were 73.3% [95% confidence interval (CI) 0.45, 0.92] and 80.0% (95% CI 0.72, 0.87), respectively, at the best cutoff point. The AUC had a value of 0.73 (95% CI 0.55–0.90) ( Figure 7C ). Multivariate survival analysis showed that the ASF1B expression level was the sole independent risk factor for the DFS rate when the ASF1B high/low expression groups were divided according to the best cutoff points determined by the ROC curve analysis [hazard ratio (HR) 7.98; 95% CI 2.54, 25.8] ( Figure 7B and Table 1 ).

A significant correlation was detected between the ASF1B expression level [fpkm log2(x+1)] and six immune factors. Spearman rank analyses demonstrated that ASF1B expression levels were significantly positively correlated with the ESTIMATE score (ρSpearman = 0.20, p < 0.05, 95% CI 0.03, 0.36), immune score (ρSpearman = 0.20, p < 0.05, 95% CI 0.03, 0.36), and T regulatory cell (Treg) infiltration level (ρSpearman = 0.22, p < 0.05, 95% CI 0.05, 0.38) ( Figure 8 ). In addition, the results based on the analysis of previously published urothelial cancer research showed that patients with high expression of ASF1B would not respond to anti-PD-L1 (atezolizumab) therapy (230 nonresponders vs. 68 responders, log2FC = 0.436, p < 0.01) (42).

The high expression of ASF1B in male TC tissues was verified in an independent male TC dataset composed of the GSE-29265, GSE-3467, and GSE-76039 datasets (independent samples t-test, p ≤ 0.01) ( Figure 9 ). These datasets were corrected for batch effects using the Combat method ( Supplementary Figure 4 ).

Three sh-RNAs were designed to silence ASF1B, and sh-ASF1B-3 (sh3) showed the best silencing effect and was used for the subsequent experiments ( Figure 10A ). CCK8 assays showed that KTC-1 cell proliferation was significantly decreased in the sh-ASF1B group than the sh-NC group, p < 0.05 ( Figure 10B ). Apoptosis assays by flow cytometry indicated that ASF1B silence increased apoptosis rate in KTC-1 cells, p < 0.05 ( Figure 10D ). Meanwhile, more cells were found to be arrested in S phase in sh-ASF1B cell lines as well, p < 0.05 ( Figure 10C ). Additionally, wound healing assay and Transwell assay were applied to evaluate male TC cell migration and invasion, respectively. We found that compared with sh-NC group, the migrated cells ( Figure 11B ) and the invasive cells ( Figure 11A ) in the sh-ASF1B group decreased significantly, p < 0.05. These results suggested that ASF1B silencing inhibited male TC cell proliferation, cell cycle progression, migration, and invasion, and induced cell apoptosis.

RT-qPCR and Western blot results showed that the silencing of ASF1B led to an increase in the mRNA and protein levels of p-AKT, AKT, and FOXP3 ( Figures 12A, B ). The results suggested that ASF1B might activate the PI3K/Akt signaling pathway through upregulating the p-Akt and Akt in male TC cells. FOXP3 also acts as a target of ASF1B, which could be involved in tumor immune escape and tumor cell proliferation.

RT-qPCR confirmed the overexpression of FOXP3 in the FOXP3-overexpressing plasmid-transfected group ( Figure 12C ; p < 0.05). In the scratch-wound assay, the width of the scratch was smaller in the FOXP3-overexpressing sh-ASF1B group than in the sh-ASF1B group ( Figures 12F, G ; p < 0.05). CCK8 assay showed that the FOXP3-overexpressing sh-ASF1B group had a significant higher proliferative ability than the sh-ASF1B group ( Figure 12E ; p < 0.05) and there was no significant difference between the FOXP3-overexpressing sh-ASF1B group and the untreated control group after a 72-h culture ( Figure 12D ; p > 0.05). These results suggest that the ability to migrate and proliferate was restored in ASF1B-silenced KTC-1 cells after being transfected with FOXP3-overexpressing plasmid.