Genome-wide association study of SNP- and gene-based approaches to identify susceptibility candidates for lupus nephritis in the Han Chinese population

Background Lupus nephritis (LN) is one of the most common and serious complications of systemic lupus erythaematosus (SLE). Genetic factors play important roles in the pathogenesis of LN and could be used to predict who might develop LN. The purpose of this study was to screen for susceptible candidates of LN across the whole genome in the Han Chinese population. Methods 592 LN patients and 453 SLE patients without renal damage were genotyped at 492,970 single nucleotide polymorphisms (SNPs) in the genome-wide association study (GWAS). Fifty-six SNPs were selected for replication in an independent cohort of 188 LN and 171 SLE without LN patients. Further quantitative real-time (qRT) PCR was carried out in 6 LN patients and 6 healthy controls. Gene-based analysis was conducted using the versatile gene-based test for GWAS. Subsequently, enrichment and pathway analyses were performed in the DAVID database. Results The GWAS analysis and the following replication research identified 9 SNPs showing suggestive correlation with LN (P<10-4). The most significant SNP was rs12606116 (18p11.32), at P=8.72×10−6. The qRT-PCR results verified the mRNA levels of LINC00470 and ADCYAP1, the closest genes to rs12606116, were significantly lower in LN patients. From the gene-based analysis, 690 genes had suggestive evidence of association (P<0.05), including LINC00470. The enrichment analysis identified the involvement of transforming growth factor beta (TGF-β) signalings in the development of LN. Lower plasma level of TGF-β1 (P<0.05) in LN patients and lower expression of transforming growth factor beta receptor 2 in lupus mice kidney (P<0.05) futher indicate the involvement of TGF-β in LN. Conclusions Our analyses identified several promising susceptibility candidates involved in LN, and further verification of these candidates was necessary.


Introduction
Lupus nephritis (LN) is a major risk factor for morbidity and mortality in systemic lupus erythaematosus (SLE), which is an autoimmune disease characterized by the loss of self-tolerance and formation of nuclear autoantibodies and immune complexes (1,2). Up to 60% of adult SLE patients will develops to LN and 10% to 30% of severe LN (Class III and above) cases might progress to end-stage renal disease (ESRD) (3). Indeed, patients with LN often have higher mortality and die earlier than the SLE patients without LN (4). LN is also associated with a 6-fold increase in mortality compared with the general population, and ESRD patients show a 26-fold increase (5). Thus, it is important to be able to predict individuals who might develop LN.
Genetic factors play an important role in the pathogenesis of LN (6). Most studies aiming to identify susceptibility genes for LN are candidate gene association analyses. In such studies, SLEassociated genes, mostly related to immunology/inflammation, are often selected to examine their role in LN. However, this approach might lead to the omission of a large number of genes with other roles in the pathogenesis of LN. Genome-wide association studies (GWAS) can explore susceptibility genes in the whole genome and have been successfully applied to multiple complex diseases.
We aimed to find susceptible candidates of LN across the whole genome in the Han Chinese population. We previously performed a GWAS for SLE by genotyping 1,047 SLE cases and 1,205 healthy controls by using Illumina Human610-Quad BeadChips in a Han Chinese population (7). Based on this GWAS data set we performed whole genome single-nucleotide polymorphism (SNP)-and gene -based analyses for the first time to explore LN susceptibility polymorphisms, genes, or loci in the Han Chinese population in the present study.

Subjects and phenotype definition
All the samples used in this study were obtained from multiple hospitals in China. Clinical information was collected from questionnaires and clinical records. Informed consent was obtained from all participants, and the study was approved by the institutional review board at each institution.
In this study, LN patients were cases, while the SLE without LN were the controls, who were genetically comparable to the LN patients. All patients met the American College of Rheumatology criteria (1997) (8). As shown in Table 1, there were no significant differences in sex or age between the two groups. LN was diagnosed according to the renal sub-phenotype data using ACR classification (proteinuria of >500 mg/24 hours or >3+, or the presence of cellular casts, or diagnosed as LN by renal biopsy). SLE without LN was diagnosed as SLE without the above renal disorder. There were 592 LN and 453 SLE without LN patients in the discovery stage, 188 LN and 171 SLE without LN patients in the replication study, 6 LN patients and 6 healthy controls in the quantitative real-time PCR experiment and 39 LN patients and 35 healthy controls in the Enzyme-linked immunosorbent assay experiment.
A total of 1,099 SLE cases were genotyped initially. Cases with genotyping call rates less than 98% or repeated or related samples were removed (52 SLE patients); 2 SLE cases without renal information were also removed. Subsequently, principal component analysis (PCA) was performed to assess genetic heterogeneity using EIGENSTRAT. No outlier was revealed by the PCA plots (S1 Figure). In total, 1,045 SLE patients (592 LN patients and 453 SLE without LN patients) were included in genome-wide association analysis.
All SNPs on the X, Y and mitochondrial chromosomes were removed, after which SNPs with unclear cluster patterns of the genotyping data were removed. SNPs with genotyping call rates less than 90% or Hardy-Weinberg equilibrium P <10 -7 were also removed. Thus, 492,970 SNPs were included in genome-wide association analysis.

Imputation
Genotypes were imputed on the basis of the 1000 Genomes Project phase 3 (hg19) using the IMPUTE program (v2.0). A total of 14,076,911 SNPs were imputed. A quality control was applied for the imputed SNPs: all SNPs on the X, Y and mitochondrial chromosomes were removed. SNPs with genotyping call rates less than 90% or Hardy-Weinberg equilibrium P <10 -7 were also removed.

SNP selection for replication analysis
As the P values of the imputed SNPs were also modest, the SNP selection for the replication study was based on the results of genome-wide association analysis without imputation. SNPs with P<5×10 -4 in the genome-wide association analysis were selected, and those with high P values of Hardy-Weinberg equilibrium (P>5×10 -4 in controls), high call rates (>95%), high minor allele frequencies (>0.05 both in cases and controls), and proximity to genes associated with LN pathogenesis (immunology, inflammation, renal resident cell, and so on) remained. The top one or two tag SNPs with the lowest P value were selected when multiple SNPs were localized to one distinct genomic locus based on physical location. Ultimately, 56 SNPs in 34 loci were selected for replication analysis.

Genotyping and quality controls in the replication study
The 56 SNPs were genotyped in 191 LN patients and 171 SLE without LN patients by using the Agena MassARRAY system. Agena MassARRAY Assay Design 3.1 software was used to design a multiplexed SNP MassEXTEND assay, and SNP genotyping was performed using Agena MassARRAY RS1000 with the manufacturer's protocols. Agena Typer 4.0 software was employed to perform data management and analyses. Three cases with call rates <90% were excluded. Two SNPs were excluded because of a call rate <90% or a P value of Hardy-Weinberg equilibrium <0.001. The remaining SNP cluster patterns from the genotyping data were checked to confirm their good quality. After quality control, 54 SNPs in 188 LN patients and 171 SLE without LN patients were remained in the replication stage analysis.

Association analysis
In the discovery stage, association analyses were performed with PLINK using logistic regression with age and sex as covariates. In the replication study, the 54 SNPs were analysed with the same association test. To combine the association evidence from the GWAS and SNP replication samples, a joint analysis of all combined samples was performed using metaanalysis with PLINK. SNPs with P<5×10 −8 were considered genome-wide significant, and those with P<10 −4 were considered as suggestive association.

Bioinformatics analysis and expression analysis
To evaluate the potential biological function of the most significant SNP, rs12606116, and related genes, we used QTLbase (http://mulinlab.org/qtlbase/), GTEx (https:// gtexportal.org/home/), and GEO (https://www.ncbi.nlm.nih. gov/geoprofiles/) databases to conduct expression quantitative trait locus (eQTL) analysis and gene expression difference analysis. The gene expression data of LN and normal kidney tissues were downloaded from the GEO database. Those genes around rs12606116 ± 400 kb were focused on. The eQTL data of rs12606116 and genes around rs12606116 ± 400 kb were downloaded from QTLbase and GTEx databases.

Gene-based and gene set enrichment analysis
Gene-based analysis was performed with the versatile genebased test for GWAS (VEGAS), which uses SNP-level data to incorporate information from a full set of markers annotated to each gene and accounts for linkage disequilibrium (LD) between markers (9, 10). The association P-values of a given gene with n SNPs were converted to upper tail chi-square statistics with 1 degree of freedom. Then the empirical gene-based P-value was calculated using the Monte Carlo simulation (11).
To explore the biological functions and cellular signaling pathways, all genes with P<0.05 in the gene-based test were used as input to perform gene set enrichment analysis using Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG). GO analysis was performed in three aspects: biological process (BP), cellular component (CC), and molecular function (MF). GO analysis and KEGG pathway analysis were both conducted by the Database for Annotation, Visualization, and Integrated Discovery (DAVID).

Animals
Three 6-8 weeks old MRL/lpr mice were purchased from Shanghai Model Organisms and subsequently maintained in the animal facility at the Chinese People's Liberation Army General Hospital. Mice were monitored for lupus for 12 weeks, at which time they were sacrificed for renal tissue. Renal tissue from three BALB/c mice were served as control. All animal experiments conducted have been approved by the Animal Ethics Committee of the Chinese People's Liberation Army General Hospital.

Quantitative real-time PCR
The gene expression differentiation of ADCYAP1 and LINC00470 genes was analysed in peripheral blood mononuclear cell (PBMC) samples of 6 LN patients and 6 healthy controls. Besides, the gene expression of transforming growth factor beta receptor 2 (TGFBR2) gene was analysed in renal tissue of 3 MRL/lpr lupus mice and 3 BALB/c mice. RNA templates were extracted from PBMC or renal tissue samples with TRIzol. cDNA was synthesized from 1 mg of RNA templates using reverse transcriptase and oligo(dT) primers (NEB). The sequences of the primers used for amplification are given in S1 Table. 18S was used as a reference for normalization. The qRT-PCR reaction was carried out in triplicate using SYBR Green on an ABI7500 Real-Time PCR System.

Enzyme-linked immunosorbent assay
Human transforming growth factor beta1 ELISA Kit (ab100647; Abcam) was used to measure plasma levels of transforming growth factor beta1 (TGF-b1) according to the manufacturer's instructions.

Statistical analysis
Data were analysed using SPSS 21.0 (SPSS, Chicago, IL, USA) statistical software. Continuous variables were expressed as median (inter-quartile range) or mean ± standard deviation according to the distribution, and categorical variables were expressed as frequency (percentage). Differences between the two groups were examined using the Mann-Whitney U-test or Student's t-test. For categorical variables, differences between groups were analysed using the c 2 test. P values<0.05 were considered statistically significant.

Genome-wide association of LN
There were 592 LN patients and 453 SLE without LN patients included in the GWAS discovery stage. As shown in Table 1, there were no significant differences in sex or age between the two groups.Manhattan and quantile-quantile plots (q-q plots) are shown in S2 Figure Table). The most significant SNP was rs10151371 (P=9.90×10 −6 ) at 14q31.3.

Imputation
A total of 14,076,911 SNPs were imputed. After quality control, the most significant imputed SNP was rs117609374 (P=3.67×10 −6 ) at 14q31.3. Considering the modest P value of the imputation results, only the imputed region plots of the 9 SNPs from 7 loci described below are presented in Figure 1.

Replication study
To assess the association of possible SNPs with LN, 56 SNPs in GWAS were selected for genotyping in another 191 LN and 171 SLE without LN patients. After quality control, 54 SNPs in 188 LN patients and 171 SLE without LN patients were analysed. Due to limited power of the replication samples, only one SNP rs12606116 exhibits P< 0.05 (S3 Table).

Meta-analysis of the GWAS and replication results
In the meta-analysis of the GWAS and SNP replication results, 9 SNPs showed a suggestive association with LN, at P< 10 -4 ( Table 2). The most significant SNP was rs12606116 (P=8.72×10 −6 ), located at 18p11.32. This SNP localizes to an intergenic region close to LINC00470, a long intergenic nonprotein coding RNA gene, and the ADCYAP1 gene, which encodes a potent renoprotective peptide. The P values of the remaining 8 SNPs were between 10 -5 and 10 -4 ( Table 2). Region plots of the 9 SNPs from 7 loci are presented in Figure 1.

qRT-PCR
Further studies were applied to verify the biological function of rs12606116 and related genes. Data in the GEO database showed that mRNA expression of ADCYAP1, which was the closest coding gene to rs12606116, was significantly different in the renal tubules of LN patients and healthy controls (GSE32591, P=1.06×10 −3 , logFC=0.16). The QTLbase database showed an association between rs12606116 and LINC00470, which was the closest gene, in CD4+ T cells and CD14+ monocytes cells. Thus, we measured the mRNA levels of LINC00470 and ADCYAP1 in PBMCs of 6 LN patients and 6 healthy controls using qRT-PCR. As shown in Figure 2, the expression of LINC00470 and ADCYAP1 in the PBMCs of LN patients was significantly lower than that in healthy controls. These results further indicated LINC00470 and ADCYAP1 in 18p11.32 might be involved in the development of LN.

Results of gene-based and gene set enrichment analysis
To explore additional susceptibility genes associated with LN, we optimally use the GWAS data sets to conducted genebased analysis using VEGAS. In the gene-based analysis, the SNPs in the GWAS were mapped to 14,012 genes (nSNP>1). In total, 690 genes showed evidence of association with LN (P<0.05), including LINC00470, identified by the above SNPbased analysis (S4 Table). Table 3 shows the top 10 genes ranked by their P-values from the VEGAS analysis. The topmost hit was LINC01146 with P=1.20×10 −5 . Two of the 10 genes also featured in the list from t h e S N P -b a s e d a n a l y s i s : CI T ( P = 2 . 78 × 1 0 − 4 ) and WDR25 (P=3.39×10 −4 ).
Furthermore, we carried out enrichment analysis with all the 690 genes those presented P<0.05 in the gene-based Region plots of the genome-wide association results of the suggestive single-nucleotide polymorphisms (SNPs analysis. The ten most significant GO terms and pathways are shown in Figure 3. Top associated GO terms included "cell adhesion", "positive regulation of MAPK cascade", "protein kinase activity", "Chemokine signaling pathway", "TNF signaling pathway", and so on. Interestingly, "transforming growth factor beta receptor signaling pathway", ''type I transforming growth factor beta receptor binding'' and ''transforming growth factor beta binding'' were significantly enriched, highlighting the involvement of transforming growth factor beta (TGF-b) signalings in the pathogenic mechanisms of LN.

Detection of related molecules in TGF-b signaling pathway
We further examined the key molecules, TGFBR2 and TGF-b1, in TGF-b signaling pathway. The MRL/lpr mouse is a spontaneous disease model for complement-associated inflammatory kidney disease, similar to LN (12,13). Thus, qRT-PCR was performed to assess the expression of TGFBR2 gene in renal tissue of 3 MRL/lpr mice (LN group) and 3 BALB/c mice (control group). As shown in Figure 4A, the expression of TGFBR2 gene in the LN group was significantly lower than that  Frontiers in Immunology frontiersin.org in controls. Then the plasma levels of TGF-b1 were detected using ELISA in 39 LN patients and 34 healthy controls. It was observed that plasma levels of TGF-b1 was significantly lower in patients with LN comparing to control group ( Figure 4B).

Discussion
We performed GWAS of SNP-and gene-based analyses to identify susceptibility SNPs, genes and loci associated with LN in the Han Chinese population for the first time. In the SNP-based GWAS results, 9 SNPs showed suggestive associations with LN (P< 10 −4 ), among which rs12606116 was the most significant SNP. Subsequently, bioinformatics analysis and qRT-PCR research verified the mRNA levels of LINC00470 and ADCYAP1, the closest genes to rs12606116, in PBMCs were significantly lower in LN patients than in healthy controls. Genebased analysis further validates the association of LINC00470 with LN. These results indicated that LINC00470 and ADCYAP1 in 18p11. 32 (14). The strongest evidence for association was observed localizing to 4q11-q13 (PDGFRA, P=4.5×10 −7 ). Chen et al. performed a meta-analysis of two GWASs of European SLE (1,152 LN and 1,949 SLE without LN patients), they also did not identify any genome-wide loci significantly associated with LN (P< 5×10 −8 ) (15). Similar to the two GWASs, we did not find SNPs with P< 5×10 −8 . The less significance of those results might be due to the control group, namely, SLE patients without renal damage. Some of these patients in the control group might already have kidney damage without clinical manifestations, or they might develop LN in the future. All these factors potentially reduce statistical power and increase type 2 error rather than false positive findings. Nonetheless, there is a clear trend in our results. We identified nine new suggestive SNPs associated with LN (P< 10 −4 ). The most significant SNP was rs12606116 (P=8.72×10 -6 ), which is located close to LINC00470 and ADCYAP1 at 18p11.32, the new locus reported in LN GWASs. Thus, we believe that the results is a new suggestive locus for LN in the Han Chinese population.
The qRT-PCR replication research in the PBMCs further illustrated our findings. In the biological analysis, we found that the mRNA expression of ADCYAP1 was different in the renal tubules of LN patients and healthy controls. The QTLbase database also showed an association between rs12606116 and LINC00470. Further qRT-PCR replicated that the expression of LINC00470 and ADCYAP1 in the PBMCs of LN patients was significantly lower than that of healthy controls. Especially, LINC00470 was also identified in the gene-based analysis.
LINC00470 is a long intergenic non-protein coding RNA gene and is affiliated with the lncRNA class. Noncoding RNAs, particularly microRNAs and lncRNAs, have been implicated in multiple biological processes in normal and various diseases. There are few reports about this gene. As reported, LINC00470 is a new AKT activator that mediates glioblastoma cell autophagy (16). It also facilitates malignant proliferation in hepatocellular carcinoma, lung cancer, and gastric cancer (17-19). However, its function in LN needs further study.
ADCYAP1 encodes pituitary adenylate cyclase-activating polypeptide (PACAP), which is a multifunctional neuropeptide that might act as a potent renoprotective peptide (20). Experiments in PACAP-deficient mice showed that a lack of endogenous PACAP leads to higher susceptibility to in vivo renal ischaemia/reperfusion, suggesting that PACAP protects the kidney against oxidative stress (21,22). In addition, PACAP may have an immunomodulatory role. As reported, PACAP was widely distributed in the organism including the immune system, such as thymus, lymph nodes and spleen. While exposed to proinflammatory stimuli, PACAP may prevent inflammation by down-regulating a series of cytokines and chemokines produced by innate immune cells (such as macrophages). It has been also shown that PACAP can participate in the immune response by inducing the production of regulatory T cells (23, 24).
Interestingly, two genes (CIT and WDR25) showed suggestive association with LN in both the SNP-based and gene-based analyses. CIT gene encodes a serine/threonineprotein kinase that functions in cell division. CIT has been reported to play a role in cancer proliferation and central nervous system development (25). WDR25 encodes a protein containing 7 WD repeats, involved in cell cycle progression, signal transduction, apoptosis, and gene regulation (26). Genebased analysis also highlighted several other genes, including a gene implicated in polycystic kidney disease (TMEM207), a gene expressed at the tight junction of renal tubule epithelial cells (CLDN1), and a gene associated with cell proliferation (GNL2) (27)(28)(29). It would be interesting to see whether these genes could be replicated in future studies.
Our enrichment analysis revealed several pathways, especially the TGF-b signaling pathways, including "transforming growth factor beta receptor signaling pathway", ''type I transforming growth factor beta receptor binding'' and ''transforming growth factor beta binding''. There have been many studies, including cohorts and animal models, suggesting the involvement of TGF-b, especially TGF-b1, in SLE. Our results also highlighted the important role of TGF-b in LN. TGF-b signalings regulate various biological responses, including proliferation, migration, differentiation, apoptosis, and the immune response (30). There are three TGF-b forms (TGF-b1, 2, and 3) and three TGF-b receptors (TGFBR1, 2, and 3). Among them, TGF-b1 was the first fully cloned of the TGF-b superfamily members, and was crucial in the signaling pathway (31). Besides, TGFBR2 is the receptor directly bounds to TGF-b, and thus it is also a key molecule for activation of downstream signaling (30). To verify the involvement of TGF-b in LN, we further detected TGFBR2 and TGF-b1. We found that TGFBR2 mRNA expression decreased in the LN mice kidneys, and plasma TGF-b1 levels also decreased in LN patients. As reported, SLE/LN patients often produce lower levels of TGF-b1 when compared with healthy individuals (32)(33)(34)(35)(36). However, there are limited studies on the mechanism of TGFBR2 in LN. Saxena et al. reported that expression of TGF-b receptor type II was increased in the kidneys of autoimmune-prone BWF1 mice (32). Differences in animal models might partly explain the opposite results. Taken together, these findings in our study indicate that reduced TGF-b production might predispose for the development of lupus and target organ damage. However, those association need to be verified in larger sample and the mechanisms remains to be further tested.
Over the past decade, an explosion of SLE GWASs in different races have greatly improved our understanding of the genetic basis of SLE. Over 100 susceptibility loci and hundreds of SNPs for polygenic, multifactorial SLE have been identified. We systematically collected 405 SLE associated SNPs from the GWAS catalogue (https://www.ebi.ac.uk/gwas/) and previous reported SLE GWASs. Among them, 144 SNPs were genotyped in this study. Results of the 144 known SLE susceptible SNPs in this LN GWAS were listed in S5 Table. Unexpectedly, only 9 SNPs showed P<0.05 in this LN GWAS, and only 3 SNPs had consistant OR (rs4622329, rs3093030, rs11117432). Furthermore, genes showed in Table 2 and Table 3 also have not been identified in SLE GWASs. It might suggest that, SNPs those play important roles in the pathogenesis of SLE might not play an equally critical role in LN. There might be some differences between SNPs involved in SLE vs. in LN. We further paid attention to the monogenic causes of SLE and the pathways involved. As reviewed by Demirkaya, there have been more than 30 genes causing the monogenic form of SLE identified, involved in pathways of complement, Type I IFN, self-tolerance, RAS, and so on (37). One patient with NS/Noonan-related syndrome phenotype and SLE characteristics had been reported to have a PTPN11 mutation (37). Coincidentally, PTPN11 gene showed a P value<0.05 in our gene-based analysis. As a member of the RAS/MAPK cascade, PTPN11 plays pivotal roles in cell proliferation, differentiation, survival and cell death (38). As we know, MAPK pathway is closely related to renal physiology and pathophysiology (39). The "positive regulation of MAPK cascade" was also enriched in our analysis (Figure 3). We speculate that PTPN11 might be involved in LN by activating MAPK downstream cytokine and pathways. The mechanism remains to be verified. Mutation in TNFAIP3, encoding the NF-kB regulatory protein A20, was also reported to associated with SLE (40). TNFAIP3 is an ubiquitinediting enzyme with a critical function in the inhibition of key proinflammatory molecules, mainly functions as an endogenous regulator of inflammation through termination of nuclear factor (NF)-kB activation (41). TNFAIP3 is downstream of TNF signaling pathway and can be induced by TNF. In our enrichment analysis, TNF signaling pathway was also enriched, indicating that TNFAIP3 and other members of TNF signaling pathway might play a role in the pathogenesis of lupus and LN (Figure 3).
There are some limitations in this study. First, the sample size in the discovery stage and SNP replication study was not appreciable enough. With 1,404 SLE participants, we had limited power, at a genome-wide significance level, to uncover an association with small-to-moderate effect sizes. Studies with larger sample sizes or multi-team cooperation of GWAS metaanalyses might be helpful for identifying significant loci. Second, as mentioned above, the control group, namely, SLE patients without renal damage, was clinically evaluated, and some of these patients might already have kidney damage without clinical manifestations. In addition, some patients might develop LN in the future. All these factors potentially reduce statistical power and increase type 2 error rather than false positive findings; thus, they do not decrease the robustness of the associations reported here. Thirdly, compared with SLE without LN patients, the heritability of LN might not be significant enough. Some variants with lower allele frequencies and weaker genetic effects that are difficult to identify might exist. Finally, the expression data of ADCYAP1 and LINC00470 genes in the kidney was not available. We have tried to analyse the gene expression differentiation of the two genes in kidney tissues of lupus prone mice. However, the expression of the two genes in kidney was too low to acquire in the qRT-PCR. Thus, we collected the PBMC, which included both the CD4+ T cells and CD14+ monocytes cells. Moreover, larger sample sizes in the qRT-PCR and ELISA parts might be needed to verify the differences.

Conclusion
In summary, we performed SNP-and gene-based GWAS of LN for the first time in the Han Chinese population and found some promising candidates that exhibit suggestive associations with LN. In particular, the association between rs12606116 in 18p11.32 and LN was replicated, indicating that it was a new suggestive locus for LN in the Han Chinese population. We also identified several other promising genes and pathways for LN, especially the involvement of TGF-b pathways. These findings enrich our understanding of LN inheritance in the Han Chinese population and provide important clues for future genetic research on LN. However, further exploration in the future is warranted.

Ethics statement
The studies involving human participants were reviewed and approved by the PLA institutional review board. Written informed consent to participate in this study was provided by the participants' legal guardian/next of kin.