A flow diagram of the study is shown in Figure 1. Briefly, we first analyzed the GSE37460 and GSE104948 and identified DEGs common to the two datasets. Hub genes were identified by examining the upregulated DEGs using the cytoHubba plugin of Cytoscape and by ROC curve analysis. Genes expressed in IgAN patients in the GSE93798 dataset were selected for WGCNA to identify key modules related to renal function (creatinine and eGFR); the associated hub genes were identified based on gene significance (GS) and module membership (MM) values. Co-hub genes with potential value for further analysis were also extracted. We then performed in vitro and in vivo experiments using mouse (M)MCs and clinical specimens, respectively, to validate the relevance of the co-hub genes to IgAN. This study was approved by the Renji Hospital Ethics Committee of Shanghai Jiaotong University School of Medicine (approval no. KY[2019]002).

The GSE37460 (Berthier et al., 2012), GSE104948 (Grayson et al., 2018), and GSE93798 (Liu et al., 2017) Series Matrix datasets were obtained by microarrays from the Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo/). The datasets were generated from renal biopsy specimens from IgAN patients and healthy control subjects (Table 1).

Differential gene expression analysis and preprocessing by normalization of the GSE37460 and GSE104948 datasets were carried out using R v4.0.0 software with the limma package (http://www.bioconductor.org/packages/release/bioc/html/limma.html). The P-value of genes in these datasets was defined as an adjusted P-value. After analysis, an adjusted P < 0.05 and |log2(fold change)| ≥1 were selected as the cutoffs for screening DEGs. A volcano plot and hierarchical clustering heatmap were used to visualize the differential expression of all genes and top 100 DEGs, respectively. DEGs present in both datasets were considered overlapping DEGs, and those that were upregulated were retained for further analysis (http://bioinformatics.psb.ugent.be/webtools/Venn/).

We used the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING, http://string-db.org/) (Szklarczyk et al., 2015) to construct a protein–protein interaction (PPI) network of the identified DEGs. The network was analyzed using the cytoHubba plugin of Cytoscape (http://apps.cytoscape.org/apps/cytohubba), which includes 12 algorithms for filtering hub genes from PPI modules. The stress and betweenness (Chin et al., 2014; Thakur et al., 2017) statistical measures were used to evaluate the importance of nodes in the PPI network. The top 15 overlapping DEGs in the integrated results obtained with the two algorithms were considered significant genes.

We used the “sva” package of R software to eliminate the batch effect on the original expression of hub genes and the “ggplot” and “ggpubr” packages to analyze the normalized expression matrix and generate a violin plot.

ROC curve analysis was performed using the “pROC” package of R software to evaluate the diagnostic accuracy of each hub gene, and the sensitivity and specificity were assessed based on area under the ROC curve (AUC). Hub genes with the top 5 AUC values were selected as key genes.

The original GSE93798 dataset included 20 IgAN patients (Supplementary Table 1), 22 healthy controls, and their clinical data and was deemed suitable for WGCNA. However, the investigators who generated the original dataset considered that only 19 samples met the research requirements of their study. We found that the age and sex of two samples derived from IgAN patients were the same, creating a mismatch between the clinical data and sample number; we therefore discarded these two samples and ultimately included 17 IgAN samples in our WGCNA.

The gene expression matrix of the 17 IgAN samples from GSE93798 was constructed, and the top 25% of genes with the largest variance in the samples were selected for WGCNA using the WGCNA package of R v4.0.0 software. Outlier samples were detected by hierarchical clustering, and an appropriate soft threshold power was selected in order to achieve standard scale-free networks. The next stage involved the construction of an adjacency matrix and topologic overlap matrix (TOM) and the calculation of corresponding dissimilarity (1 – TOM); a gene dendrogram was generated and module identification was performed by dynamic tree cutting with the minimum module size set to 30. The clustering of module eigengenes was carried out to merge highly similar modules with a dissimilarity of 0.25. The gene network was constructed with 1,000 randomly selected genes.

Module eigengenes were defined as the first principal component of the expression matrix for a given module and represented the average gene expression level for all genes in each module. Demographic and clinical information including sex, age, disease stage, and renal function were converted to numerical values, and regression analysis was performed between these and the module eigengene values. MM was defined as the association between a gene and a given module, and GS was defined as the correlation between genes and clinical features. Genes with high GS and MM for a clinical feature were retained as candidates for subsequent analysis. All analyses were performed using the WGCNA package of R v4.0.0 software.

We used the Database for Annotation, Visualization and Integrated Discovery (DAVID; http://david.abcc.ncifcrf.gov/) (Huang da et al., 2009) to perform Gene Ontology (GO) function and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses for overlapping DEGs, with P < 0.05 and gene count ≥2 used as cutoffs. The modules of interest were visualized using STRING (string-db.org/cgi/input.pl).

Antibodies against periostin (POSTN; sc-398631) and glyceraldehyde 3-phosphate dehydrogenase (sc-66163) were purchased from Santa Cruz Biotechnology (Santa Cruz, CA, USA). Recombinant mouse transforming growth factor β1 (TGF-β1) (P04202) was from MedChemExpress (Princeton, NJ, USA). Antibodies against proliferating cell nuclear antigen (PCNA; #13110) and B cell lymphoma-associated X protein (BAX; #2772) were from Cell Signaling Technology (Danvers, MA, USA). Recombinant mouse POSTN/osteoblast-specific factor (OSF)-2 protein (2955-F2-050) was from R&D Systems (Minneapolis, MN, USA). Cell Counting Kit 8 (CCK-8) was from Dojindo (Shanghai, China). The terminal deoxynucleotidyl transferase dUTP nick end labeling (TUNEL) assay kit was from Roche (Basel, Switzerland). Horseradish peroxidase-conjugated goat anti-mouse IgG H&L (ab6789) and goat anti-rabbit IgG H&L (ab205718) were from Abcam (Cambridge, MA, USA). Small interfering (si)RNA against POSTN and a scrambled negative control (NC) siRNA were obtained from Genomeditech (Shanghai, China) and had the following sequences: POSTN, 5′-CAAUGUCAAUGUUGAGCUA UAGCUCAACAUUGACAUUG-3′ and scrambled, 5′-UUCUCCGAACGUGUCACGUdTdT ACGUGACACGUUCGGAGAAdTdT-3′. POSTN RNA primers were synthesized by Sangon Biotech (Shanghai, China) and had the following sequences: forward, 5′-AAGGGAATGACTAGCGAAGAAA-3′ and reverse, 5′-CTCATTCACTAGAAGCGTTTCG-3′.

MMCs were obtained from the Cell Bank of the Chinese Academy of Sciences (Shanghai, China) and cultured in Dulbecco's modified Eagle's medium-F12 with 10% fetal bovine serum. After serum starvation for 24 h, the cells were stimulated with different concentrations (1, 2.5, 5, 10, and 20 ng/ml) of TGF-β1 for 48 h and recombinant mouse POSTN protein for 24 h. siRNA against POSTN was used to knock down POSTN gene expression in MMCs, which were then stimulated with TGF-β1 (20 ng/ml) for 48 h.

Western blotting and real-time (RT)-PCR were used to detect protein and mRNA levels of POSTN, respectively. BAX and PCNA levels were also evaluated by western blotting. Cell proliferation was assessed using CCK-8, and apoptotic cells were detected with the TUNEL assay.

We enrolled 120 inpatients who underwent renal biopsy at the Nephrology Department of our hospital from June 2019 to June 2020 after the patients provided informed consent to participate in the study. We collected clinical data and renal biopsy tissue specimens from 41 patients diagnosed with IgAN. The exclusion criteria were as follows: patients with diabetes, hepatitis, cirrhosis, severe metabolic syndrome, systemic lupus erythematosus, and secondary IgAN. According to these criteria, five patients with hyperlipidemia and one with hepatitis B were excluded; therefore, 35 patients were ultimately included in the study (Supplementary Table 2). Additionally, 15 patients (Supplementary Table 3) with kidney cancer were recruited at the Department of Urology of our hospital, and kidney tissue samples were obtained from a site 5 cm beside the cancer in these patients and used as control tissue.

Immunohistochemistry (IHC) was performed to analyze POSTN and PCNA expression in paraffin-embedded renal tissue sections. TUNEL was used to detect apoptotic glomerular mesangial cells. Two pathologists blinded to group assignment counted cells in five randomly selected fields in areas of positive staining. A semiquantitative analysis was performed using Image-Pro Plus v6.0 software (Media Cybernetics, Rockville, MD, USA) to assess protein levels of hub genes in the glomerulus region based on average optical density measurements (integrated optical density/glomeruli area).

Prism 7 software (GraphPad, La Jolla, CA, USA) was used for data analysis. Data are expressed as mean ± SEM. Differences between means were evaluated by analysis of variance and the two-tailed unpaired Student's t-test. P < 0.05 was considered statistically significant.

After normalization of GSE37460 and GSE104948 datasets (Supplementary Figure 1), we performed a screen of DEGs using the R software. In the GSE37460 dataset, we identified 181 DEGs based on the cutoff values, including 104 genes that were upregulated and 77 that were downregulated in IgAN samples relative to healthy controls (Supplementary Table 4). In the GSE104948 dataset, there were 199 DEGs, including 138 that were upregulated and 61 that were downregulated compared with the control (Supplementary Table 5). All DEGs and hierarchical cluster (heatmap) analysis of the top 100 DEGs in the two datasets are shown in Figures 2A–D. There were 37 upregulated and 21 downregulated DEGs that overlapped between the GSE37460 and GSE104948 datasets (Figures 2E,F and Table 2). GO and KEGG analyses of the 37 upregulated overlapping DEGs are shown in Supplementary Figure 2 and Supplementary Tables 6, 7.

The top 15 hub genes and sub-PPI modules were identified using cytoHubba (Figures 3A,B). We compared the hub genes and the top 15 DEGs identified using the stress and betweenness algorithms of cytoHubba (Supplementary Table 8) and found that they matched.

The results of the original gene expression analysis of the 15 hub genes were visualized using a violin plot (Figure 3C and Supplementary Table 9). We performed an ROC curve analysis and compared the AUC values of the 15 hub genes to evaluate their sensitivity and specificity for diagnosing IgAN. Five of the hub genes [POSTN, GATA-binding protein (GATA3), complement C1q A chain (C1QA), acetyl-coenzyme A acetyltransferase 2 (ACAT), and fibrinogen-like protein (FGL2)] had an AUC >0.90 and were retained as the genes with the highest potential diagnostic value for IgAN (Figure 3D and Table 3). To confirm their clinical utility, we performed an ROC curve analysis with the GSE37460 dataset; POSTN showed the highest AUC (91.8%; 95% CI: 82.8–100.0%) (Figure 3E).

WGCNA is a method for analyzing the expression patterns of multiple genes in different samples; genes with similar expression patterns are grouped into clusters or modules (Beckerman et al., 2017) and likely have the same biological functions. This allows the study of associations between modules and sample characteristics such as clinical features (Bao et al., 2020).

A total of 17 samples with clinical data were included in the WGCNA (Figure 4A). We found 4,926 genes with altered expression in IgAN relative to control samples, and the top 25% were retained for further analysis. We selected β = 7 (scale-free R² = 0.95) as the soft-thresholding power to obtain a scale-free network (Figures 4B,C). We identified 30 gene coexpression modules after excluding the gray module using a merged dynamic tree cut (Figure 4D). The heatmap of the TOM of 1,000 selected genes revealed that each module independently validated the others (Figure 4E). We also found that the module containing POSTN (i.e., the royalblue module) had the highest correlation with creatinine (R² = 0.54, P = 0.02) and eGFR (R² = −0.45, P = 0.07) (Figure 4F).

We generated scatterplots of GS vs. MM in the royalblue module with clinical information of IgAN patients. The GS values for creatinine and eGFR showed correlations of 0.53 (P = 1.2e−08) and 0.21 (P = 0.035), respectively (Figures 5A,B). Meanwhile, POSTN had very high MM and GS for creatinine (MM = 0.82 and GS = 0.66) and eGFR (MM = 0.82 and GS = −0.67) (Supplementary Table 10). Thus, POSTN was deemed as a key gene in the royalblue module. A heatmap and bar graph of gene expression distribution in the royalblue module are shown in Figure 5C.

We used the DAVID online database to perform GO and KEGG enrichment analyses of genes in the royalblue coexpression module. The genes were significantly enriched in the following biological process GO terms: cell adhesion, extracellular matrix organization, response to lipopolysaccharide, cellular response to transforming growth factor beta stimulus, and others with the threshold set to P < 0.05 and count ≥2 (Figure 5D and Table 4). There were four enrichment KEGG pathways (protein digestion and absorption, ECM–receptor interaction, focal adhesion, and amoebiasis) of these genes (Figure 5E and Table 5).

The PPI analysis of the genes in this module revealed POSTN at the core of the PPI network (Figure 6A). The stress (Figure 6B) and betweenness (Figure 6C) algorithms in cytoHubba confirmed that POSTN had high connectivity with other genes in the module.

To validate the results of the WGCNA, we carried out in vitro experiments using MMCs. Exposure to different concentrations (1, 2.5, 5, 10, and 20 ng/ml) of TGF-β1 for 48 h significantly increased the level of POSTN (Figures 7A,B), with the highest expression observed at 20 ng/ml. Therefore, this concentration was used to stimulate MMCs in subsequent experiments.

Stimulation with TGF-β1 (20 ng/ml) for 48 h significantly increased POSTN and PCNA and decreased BAX expression (Figures 7C–G). Treatment of MMCs with different concentrations (1, 2.5, 5, 10, and 20 ng/ml) of recombinant POSTN protein for 24 h had the same effects on PCNA and BAX levels as TGF-β1 (Figures 7H,I). The results of the CCK-8 assay showed that cell proliferation was increased in a concentration-dependent manner upon exposure to recombinant POSTN; the rate of proliferation was the highest in cells treated with 20 ng POSTN (Figure 7J), whereas at this concentration, the rate of apoptosis was decreased (Figure 7K).

We used POSTN siRNA to silence POSTN in MMCs for 24 h, followed by stimulation with TGF-β1 (20 ng/ml) for 48 h. POSTN depletion reduced POSTN and PCNA levels and increased that of BAX (Figures 7L–P). Cell proliferation was also decreased upon siRNA-mediated knockdown of POSTN (Figure 7Q), whereas apoptosis was increased relative to cells transfected with NC siRNA and treated with TGF-β1 (Figure 7R).

We validated the clinical relevance of our findings in MMCs by IHC using renal tissue samples from IgAN patients and healthy control subjects (Table 6). POSTN and PCNA levels were increased, whereas the rate of apoptosis detected by TUNEL was decreased in IgAN patient samples compared with those of healthy subjects (Figures 8A–D). Moreover, POSTN level in the renal tissues of IgAN patients was positively correlated with serum creatinine level (r = 0.8603, R² = 0.7401, 95% CI: 0.7388–0.9276) and negatively correlated with eGFR [r = −0.8673, R² = 0.7522, 95% CI: (−0.9314)–(−0.7511)], suggesting that renal dysfunction in IgAN involves POSTN (Figures 8E,F).