Circulating miRNAs in Type 2 Diabetic Patients with and without Albuminuria in Malaysia

Diabetic kidney disease (DKD) accounts for 39% of end-stage renal disease occurrences in the USA [1]. Likewise, DKD continues to be the leading cause of chronic kidney disease (CKD) in Malaysia, contributing to 69% of new dialysis patients in 2018 [2]. The pathology of DKD is thought to be driven mainly by hyperglycaemia, resulting in activation of intracellular signalling pathways and ultimately renal fibrosis [3]. Despite advances in management, current biomarkers and therapies fail to prevent its progression to end-stage renal disease [3]. Hence, a better understanding of the mediators and mechanisms of DKD is urgently needed.

Microribonucleic acid (miRNA) is a family of small non-coding RNA which has emerged as important post-transcriptional regulators of gene expression [4]. Aberrant expression of miRNAs and their involvement in cellular mechanism of DKD have been observed [3]. This study aimed to identify the pattern of expression of selected miRNAs in the sera of type 2 diabetes mellitus (T2DM) patients with and without albuminuria. Worldwide, similar studies have been done in various populations, as previously reviewed [5-7]. We aimed to investigate the pattern of differential expression of a selection of miRNAs which are associated with fibrotic pathway in T2DM patients, some of which have been previously studied in DKD. Previous studies have suggested ethnic-related differences in miRNA expression [8, 9]. Therefore, in addition to being able to confirm previous findings, this study may potentially reveal novel miRNAs, perhaps uniquely expressed in the Malaysian population.

This cross-sectional study involved 41 existing T2DM patients on follow-up at a community clinic, from January to April 2019. This study involved a 2-step procedure: (1) miRNAs in the sera of 12 patients (screening cohort) were profiled and (2) the most relevant miRNAs were validated with the sera from 33 patients (referred to as the validation cohort). This included 4 patients from the screening cohort. Written consent was obtained from individual patients. Ethical approval for this study was obtained from the Medical Research and Ethics Committee (MREC), Ministry of Health Malaysia. Inclusion criteria were T2DM patients and available clinical and laboratory information. Patients with type 1 DM; non-DKDs; comorbidities such as cancer, acute febrile illness, or recent infections; pregnancy; and postpartum were excluded.

Information on the following were retrieved from patients’ medical records: (1) demographic data including age and gender; (2) clinical data such as blood pressure, weight, duration of T2DM, treatment, and other diabetic complications; (3) laboratory data including fasting plasma glucose, glycated haemoglobin A1c (HbA1c), albumin-to-creatinine ratio, and serum creatinine. Patients were divided into 3 groups according to their albuminuria levels, normoalbuminuria (NA), microalbuminuria (MIC), and macroalbuminuria (MAC).

In this study, NA, MIC, and MAC were defined as an albumin-to-creatinine ratio of <30, between 30 and 300 and >300 mg/g creatinine, respectively, in at least 2 of 3 spot urine collections [10]. The estimated glomerular filtration rate (eGFR) was calculated based on the Chronic Kidney Disease-Epidemiology creatinine equation.

Collected samples of fasting (≤8 h) whole blood in 10-mL BD Vacutainer®SST^TM tubes (Becton Dickinson, Franklin Lakes, NJ, USA) were left at room temperature for 30-min clotting before 10-min centrifugation at 1,900 g. The upper phase was harvested, filtered with 0.22-μm Minisart®PES filter (Santorius, Gottingen, Germany), and aliquoted into microcentrifuge tubes for storage at −80°C until analysis.

Total miRNA was isolated using the miRNeasy Serum/Plasma Kit (Qiagen GmbH) according to the manufacturer’s protocol. In brief, the Qiazol Lysis Reagent was added to an aliquot of 200 μL of serum for denaturation followed by miRNeasy Serum/Plasma Spike-In Control C. elegans miR-39 miRNA mimic (Qiagen GmbH) to monitor RNA recovery and reverse transcription efficiency. Chloroform was added for phase separation followed by 100% ethanol addition to the upper aqueous phase. Rinse steps were achieved by Buffer RWT, Buffer RPE, and 80% ethanol. RNA was eluted by adding 14 μL RNase-free water to the RNeasy MinElute spin column membrane. RNA quantity and quality were determined by a NanoDrop 1000 spectrophotometer (Thermo Fisher Scientific, Inc., Wilmington, DE, USA), before storage at −80°C until used. Reverse transcription (RT) of miRNA to cDNA and preamplification were performed with the miScript RTII kit (Qiagen GmbH) and miScript PreAMP PCR Kit (Qiagen GmbH), respectively, according to the manufacturer’s protocol, using 6 μL template RNA and 5 μL of diluted cDNA, respectively.

Circulating miRNAs in the serum of the screening cohort were identified by the pathway-focused (human fibrosis) miScript miRNA qPCR array (Qiagen GmbH), and selected miRNAs were later validated in the validation cohort with a custom miScript miRNA qPCR array (CMIHS02742) (Qiagen GmbH) according to manufacturer’s instruction. Each reaction included cDNA template, 2x QuantiTect SYBR Green PCR Master Mix, 10x miScript Universal Primer, 10x miScript Primer Assay and RNase-free water, and diluted template cDNA to make up a total volume of 25 μL. Quantitative real-time PCR was performed on Eppendorf Mastercycler ep realplex 4, which was programmed as follows: 15 min at 95°C, 15 s at 94°C, 30 s at 55°C, and 30 s at 70°C with 40 cycles.

Real-time PCR data analysis was performed using Qiagen GeneGlobe Data Analysis Center software (https://geneglobe.qiagen.com/analyze/). Relative expression (RE) was derived by the 2^−∆Ct method, whereby expression of each miRNA was normalized to that of the housekeeping gene SNORD95 or Global Ct Mean. Fold change was calculated from the normalized miRNA expression in each group (MIC or MAC) divided by the normalized miRNA expression in the control group (NA). p values were calculated based on an independent samples t test of the replicate RE (2^−∆Ct) values for each gene in the control group and test groups. Further verification of qPCR data analysis was done manually using Microsoft Excel and IBM SPSS (Statistical Package for the Social Science; IBM Corp, Armonk, NY, USA) release 25 for Microsoft Windows software.

Initially, relative quantification of miRNAs in the screening cohort was based on the global Ct mean normalization method, and Ct > 30 was considered as negative expression, as recommended by the Qiagen GeneGlobe Data Analysis Center (https://geneglobe.qiagen.com/analyze/). Online supplementary Table 1 (for all online suppl. material, see www.karger.com/doi/10.1159/000518866) shows a RE of 84 miRNAs in the screening cohort normalized to the global Ct mean. MiRNAs with significant ±2.0-fold changes and of good quality were chosen for validation, as defined by the software. In addition, other miRNAs were chosen based on the ranking of their p values and their involvement in previous studies in DKD and DKD-related pathways, whereby miRNAs with the highest number of putative messenger RNA targets were chosen. Online supplementary Table 2 shows the additional criteria used for selection of miRNAs for validation. In validation analysis, miRNA expression was normalized to SNORD95 as the global Ct mean method was not suitable and Ct > 35 was considered as negative expression. To standardize, we re-analysed the RE of the selected miRNA in the screening cohort based on SNORD95 normalization, and this was used for subsequent analysis in comparison with the validation cohort.

Normality check of continuous variables was done using both statistical tests (Shapiro-Wilk, Kolmogrov-Smimov, skewness, and kurtosis) and visual inspection of plots (histogram, normal Q-Q, and detrended normal Q-Q). Continuous data were presented as either mean ± standard deviation or median (interquartile range). Comparison of 2 independent groups was done using independent samples t test (Tables 2–4). Comparison of >2 independent groups was done using the independent samples Kruskal-Wallis test (Table 1). Comparison of categorical variables was done using either χ² test of contingencies or Fisher’s exact test (Table 1). Correlation between selected miRNA RE and continuous variables (demographical and clinical data) was done using Pearson’s correlation coefficient. Multiple Regression test was used to analyse the linear relationship between selected miRNA RE and its significant predictors. Statistical calculations were done using IBM SPSS (Statistical Package for the Social Science; IBM Corp) release 25 for Microsoft Windows and the statistical software programme package Prism 8 (GraphPad, San Diego, CA, USA). The statistical significance was defined as p < 0.05.

There was no statistically significant difference of demographic and clinical parameters between the 3 groups of the study population, except in the levels of eGFR, urine protein, total cholesterol (TC), and usage of oral-hypoglycaemic agents and insulin (Table 1). Post hoc analyses showed that T2DM patients with MAC had significantly lower eGFR (p = 0.01) than those with NA. Significantly higher urine protein was also seen in MAC (p < 0.001) and in MIC (p < 0.002) than those with NA.

During data analysis in the initial screening cohort normalized to the global Ct mean, 12 serum samples were included (NA = 5, MIC = 3 and MAC = 4). In the reanalysis of the screening cohort normalized to SNORD95, 3 samples were excluded (1 sample from the NA group due to SNORD95 Ct value >35; 1 sample from the NA group due to poor PCR array reproducibility, and 1 sample from the MAC group due to poor RT efficiency). In the validation cohort normalized to SNORD95, 2 samples were excluded from the MIC group and 1 sample from the MAC group due to a SNORD95 Ct value >35.

Tables 2-4 show qPCR data analysis obtained from Qiagen GeneGlobe Data Analysis Center software. Results in Tables 3 and 4 tallied with those from manual verification. However, those in Table 2 did not, despite using the methods used by the software. We reported the results by the software for the following reasons: (1) the software is developed by the manufacturer; therefore, it will be the best suited and optimized for their products; (2) verification of the qPCR data normalized to SNORD95 by manual calculation supports that the software was correct in their analysis; (3) there are certain aspects of qPCR data analysis which we can only obtain from the software and not from manual analysis, whereby miRNAs with significant ±2.0-fold changes and of good quality were chosen for validation, as defined by software; (4) we have encountered problems with normalization method using the global Ct mean and have therefore reanalysed our screening cohort using SNORD95 normalization and used this for comparison with validation stage and conclusion of our study. We have however included Table 2 and Figure 1A (I)–C (I) to show that our initial selection of differentially expressed miRNAs for validation were based on the global Ct mean normalization method, as recommended by the software.

In the initial screening cohort normalized to the global Ct mean, 17 out of 84 miRNAs in the panel were chosen for validation study in serum samples. Three differentially expressed miRNAs showed statistical significance, namely, miR-142-3p (p < 0.001), miR-34a-5p (p = 0.01) in MIC, and miR-874-3p (p = 0.04) in MAC compared to NA (Table 2).

Following re-normalization with SNORD95 in the screening cohort, an upregulation of miR-142-3p in the MIC group and that of miR-874-3p in the MAC group change did not reach statistical significance, with p = 0.09 and p = 0.11, respectively. However, an upregulation of miR-34a-5p in the MIC group was statistically significant (p = 0.03). In addition, there was a statistically significant upregulation of miR-142-3p in the MAC group (p = 0.02). Statistically non-significant upregulation of miR-145-5p and miR-101-3p in the MAC groups was also observed (Table 3).

Following the validation study, the upregulation of miR-142-3p and miR-34a-5p in MAC and MIC groups, respectively, were <2-fold changes, of poor quality and statistically non-significant. Meanwhile, miR-874-3p, miR-101-3p, miR-744-5p, and miR-661 showed statistically significant differential expression of >2-fold (Table 4).

Figure 1 shows RE of (A) miR-874-3p; (B) miR-101-3p; and (C) miR-145-5p in (I) the screening cohort normalized to the global CT mean; (II) the screening cohort normalized to SNORD95; and (III) the validation cohort normalized to SNORD95. Although upregulations of miR-874-3p and miR-101-3p in the screening cohort normalized to SNORD95 (Table 3) were statistically non-significant, they were significant in the validation cohort (Table 4). Furthermore, significant upregulation of miR-874-3p was initially observed in the screening cohort normalized to the global Ct mean (Table 2). Meanwhile, miR-145-3p showed trends of upregulation in both cohorts normalized to SNORD95, albeit of statistical non-significance (Tables 3, 3, 4). For these reasons, we have concluded that a trend of upregulation for these miRNAs was seen in serum of T2DM patients with MAC compared to NA.

Pearson’s correlation coefficient analysis indicated that there was a significant medium negative correlation with the eGFR (p = 0.05) and a medium positive correlation between miR-874-3p RE and age (p = 0.03) (shown in Fig. 2). Similarly, significant medium negative correlations were seen between the eGFR and miR-101-3p (shown in Fig. 3) and miR-145-5p (shown in Fig. 4), with p = 0.03 and p = 0.05, respectively.

In this study, in keeping with the underlying renal fibrosis in DKD, we had first screened 84 miRNAs involved in human fibrosis pathways. The findings of statistical significance only in validation and not in screening stage led us to conclude that trends of upregulation of several miRNAs were observed in this preliminary study. These included upregulation of previously known miRNA, namely, miR-145-5p and possibly novel ones, namely, miR-874-3p and miR-101-3p in the serum of T2DM patients with MAC compared to those with NA.

The study population comprised confirmed T2DM patients who were on regular follow-up at a community clinic. Majority of the study population consisted of Malay ethnicity, whereby only 1 non-Malay (Chinese) patient was included. The Malay makes up the largest ethnic group [11] with the second highest prevalence of T2DM amongst the 3 main Malaysian ethnic groups of Malay, Chinese, and Indians [10]. Grouping into albuminuria groups were strictly based on 2 out of 3 positive results of their urine albuminuria levels, according to the definition implemented in clinical setting [10]. Only patients with persistent MIC levels were selected and those with intermittent MIC levels were not. This was to minimize any bias in interpretation of results based on albuminuria levels.

The study population was comparable in terms of demographical characteristics, including age and duration of T2DM and clinical parameters. However, there was statistically significant difference in the prescription of glucose-lowering medication and insulin across albuminuria groups, which was also reported in previous urine study [12]. This may reflect the increasing need for a tighter glycaemic control by means of insulin, in order to achieve a similar degree of glycaemia as patient progress from NA to MAC. Such a similar degree of glycaemia across the groups was reflected by statistically non-significant differences in fasting plasma glucose and Hba1c. There was no significant difference in most of clinical characteristics except in eGFR levels, urine protein levels, and TC levels. Previous studies have shown similar significant differences in serum creatinine and eGFR among different groups of albuminuria [13-17].

In our initial screening cohort (global Ct mean normalization), miR-142-3p, miR-34a-5p, and miR-874-3p showed significant upregulation. Other selected miRNAs included some well-known miRNAs previously studied in the sera, such as miR-29a-3p, miR-29b-5p, miR-21-5p, and miR-126-3p and a few were potentially novel, such as miR378a-3p and miR-661. Several of the well-known DKD-related miRNAs, such as miR-192 and miR-126, showed contradictory expression patterns in our study. In the initial screening, we found no dysregulation of miR-192-5p; thus, it was not selected for validation. Mixed patterns of dysregulation of circulating miR-192 of upregulation [15, 18, 19] and downregulation [16, 20-22] in MAC compared to NA have been previously shown. Similar mixed dysregulation patterns, of upregulation [6, 23] and downregulation [6, 14], were also shown for miR-126-3p. mir-126-5p was upregulated in our initial screening, albeit of statistical non-significance but subsequently showed no dysregulation in the validation cohort. Several selected miRNAs showed similar dysregulation patterns in the initial screening cohort, as previously showed, albeit of statistically non-significance. For instance, upregulations of miR-21-5p and miR-29a-3p concurred with previous findings [15, 24].

We found no dysregulation patterns of other glomeruli-enriched miRNAs previously shown in experimental models of kidney injury, namely, miR-200 family, -217, -216, and -377, as previously reviewed [25], in our initial screening cohort, and therefore, these were not selected for validation. In another review, downregulation of glomeruli-related miRNAs, namely, miR-25, -26a, -30b, -141-3p, -196a, and -let-7 and upregulation of miR-744, -146a, and -34a were shown in experimental models of DKD, as previously reviewed [26]. In line with one of these findings, miR-25 has been shown to be downregulated in the sera of T2DM patients with and without DKD [27]. In our initial screening analysis, a statistically non-significant upregulation of miR-25-3p was seen in MIC and MAC groups compared to NA. Also, similar differential expression patterns of miR-141-3p, 196a-5, and -744, were seen, although they were of statistically non-significance except for miR-34a-5p. In contrast, our initial screening analysis showed an upregulation of miR-26a-5p.

These not uncommonly reported opposing findings may be accounted for by differences in inclusion criteria, sample selection, and methodological factors such as differences in RNA extraction, miRNA quantification, and normalization methods used. Nevertheless, the difference observed could be attributable to the heterogeneity of the study population, possibly reflecting ethnic-specific differences in miRNA expression, as previously suggested [8, 9]. Supporting this notion is the observation that South Asians have unique metabolic profile which may predispose them to higher risk of T2DM than Caucasians [8]. This was further supported by a previous study which showed that the association between plasma miR-144 expression and T2DM differed between Swede and Iraqi ethnic groups [9]. Furthermore, the lack of conclusive genetic candidate accountable for ethnic variation and familial clustering of DKD may point towards epigenetic candidacy instead [28].

We concluded that a trend of upregulation of miR-874-3p, miR-101-3p, and miR-145-5p was observed in this study.

Previous studies on expression levels of miR-874-3p have shown mixed results. To the best of our knowledge, dysregulation of miR-874-3p has not been documented in any studies on circulating miRNAs in DKD. Downregulation of the circulating plasma level of miR-874-3p was however shown in obese and overweight children, associated mainly with thyroid hormone signalling and tyrosine metabolism pathways and was thought to be a potentially useful biomarker of pre-DM [29]. The upregulation of miR-874-3p in this study concurred with in vitro studies in renal podocytes and in vivo mice models of DKD, whereby miR-874-3p exacerbated doxorubicin-induced renal podocytes injury via apoptosis and oxidative damage by targeting methionine sulfoxide reductase B3 [30]. In a related renal interstitial fibrosis model, upregulation of miR-874-3p was seen in the sera of mice following aristolochic acid-induced renal injury, implicating mitogen-activated protein kinase signalling pathway [31]. Interestingly, an upregulation of miR-874-3p was seen in the kidneys of non-diabetic mice [32], although its clinical significance was not discussed by the authors.

However, in contrast to our finding, downregulation of miR-874-3p was observed in the in vitro mouse podocytes model, leading to upregulation of its target toll-like receptor 4 resulting in overexpression of inflammatory cytokines such as interleukin-, tumour necrosis factor-α, and IL-1b [33]. Similar downregulation of miR-874-3p was shown in hyperglycaemia-treated human kidney-2 cells, which was reversed by administration of psoralen, leading to reduction in inflammatory cytokines and extracellular accumulation toll-like receptor 4/nuclear factor kappa beta and Smad2 signalling pathway [34]. Downregulation of miR-874-3p was also reported in non-DKD conditions, such as in oesophageal squamous cell carcinoma tissues and cell lines, implicating signal transducer and activator of transcription (STAT)-3 as its downstream target, affecting cell proliferation, migration, and invasion [35].

Similarly, to the best of our knowledge, an upregulation of circulating miR-101-3p in patients with DKD has not been shown. In studies related to DM however miR-101 was previously shown to be significantly increased in in the sera of T2DM patients compared to non-diabetic patients [36]. The study showed that miR-101 was positively correlated with age, body mass index, HbA1c, and postprandial glucose levels, whilst HbA1c was an independent predictor of miR-101 [36]. However, correlation between miR-101 and the eGFR was not done. A prospective study showed that miR-101 could predict T2DM remission in patients with coronary heart disease [37]. However, the authors did not discuss the pattern of its dysregulation in this study [37]. An upregulation of circulating miR-101-3p was observed in T1DM patients [38, 39], and its upregulation was thought to precede the disease onset [39]. In vivo upregulation of miR-101-3p was seen in umbilical vein endothelial cells derived from patients with gestational diabetes mellitus [40], whilst upregulation of circulating miR-101-3p was also reported in non-diabetic patients such as acute coronary syndrome [41]. In contrast, an in vitro mice model of renal fibrosis showed downregulation of miR-101a in kidney tissues and human renal tubular epithelial cells [42].

With regards to miR-145-5p, our findings concurred with several in vivo studies in DKD and CKD patients, involving the sera and plasma as follow. Interestingly, although it was done in a cohort of CKD patients, with an unknown number of DKD patients, a similar statistically non-significant upregulation of miR-145 in the sera of these patients had been previously reported [43]. An increased expression of miR-145 in the plasma from T1DM patients DKD compared with healthy control subjects has also been previously shown [44].

Although the following studies were done in the urine, interestingly, upregulations of miR-145 were shown, in keeping with our findings in the sera. An in vivo upregulation of miR-145 in urinary exosomes from T1DM patients with MIC and in the mice model of early DKD has been observed [45]. Similar upregulation was seen ex vivo in mice glomeruli and in vitro in mesangial cells and their associated exosomes [45]. Despite measuring urinary total miRNA, similar upregulation of miR-145-5p was also shown in T1DM patients with MIC, which interestingly, preceded MIC by several years [46]. Subsequent study also showed an upregulation of total urinary miR-145 in CKD patients with fast compared to slow progression of the disease [47]. In contrast, although a downregulation of miR-145 in serum of stage III-IV CKD patients has been observed, the authors were analysing vascular remodelling in the context of CKD [48].

Of interest, it has been previously suggested that opposite modulation of each miRNA observed in different cells, as discussed here, may reflect varying miRNA modulation in different cell types [49]. Also, this observation is in keeping with the inherent ability of each miRNA to target several messenger RNAs [4].

The significant negative correlations between miR-874-3p, miR-101-3p, and miR-145-5p with the eGFR indicated that their patterns of differential expression may contribute towards impairment of renal function in these patients. Stronger correlation between a differentially expressed miRNA in DKD patients with the eGFR had been previously shown [50]. In contrast to previous studies [13, 20] however, in this study, no correlation was seen between levels of urinary protein with any of the differentially expressed miRNAs. This could be due to the semi-quantitative method of measuring urine proteinuria and MIC.

Unfortunately, we have encountered an issue of normalization method in our study. Because the global Ct mean normalization method was not suitable in the validation cohort as it involved using one of the miRNA as a normalizer, we re-analysed the selected 17 miRNA expression in the screening cohort to the SNORD95 normalization method. We have chosen a higher Ct value for normalization with SNORD95 compared to that used for the global Ct mean as to include as many samples as possible during data analysis. Differences observed between the 2 normalization methods in this study shows that the choice of normalization method used may partly be responsible for the incongruent results seen in different miRNA studies. At the moment, the lack of standardized normalization method is an acknowledged limitation of miRNA expression studies [51].

This study was limited by the relatively small sample size. To the best of our knowledge, this was the first miRNA profiling study in Malaysia. To overcome inherent biological heterogeneity, ideally, a repeat study involving a larger sample size should be done. Furthermore, it would be interesting to repeat this in other ethnic groups in Malaysia as this study involved mainly the Malay ethnic group.

In this study, efforts were taken to ensure that standardization of sample collection, preparation, and downstream experiments were preserved. However, a future study involving better sample preparation, for example, extracting exosome and using a different RNA extraction method may be indicated. Exosome, although more laborious, has been shown to carry purer miRNA [52].

Certain inherent issues of miRNA quantifications were encountered in this study, such as the lack of standardized normalization method. In addition, differences in platforms, study models, and samples utilized in previous studies may account for some discordance from ours and amongst them. Nevertheless, our findings of differentially expressed miR-874-3p and miR-101-3p may be unique in the Malaysian population, perhaps highlighting the role of ethnicity as previously suggested [8, 9]. In conclusion, the identification of both previously known and potentially novel differentially expressed miRNAs in DKD patients in this study has provided encouraging evidence to further investigate the putative roles of these differentially expressed miRNAs.

The authors would like to thank the Director General of Health Malaysia for the permission to publish this paper. We also thank the staff of the community clinic at Bandar Seri Putra for their support in this study. This study was funded by the UPM Grant (GP/2018/9641100) and UPM Postgraduate Grant (GP-IPS/2019/9678300).

Written informed consent was obtained from individual patients. Ethical approval was obtained from the MREC, Ministry of Health Malaysia (NMRR-17-2715-38571 [IIR]).

The authors have no conflicts of interest to declare.

This study was funded by the UPM Grant (GP/2018/9641100) and UPM Postgraduate Grant (GP-IPS/2019/9678300).

S.Y.Z.S. and H.F.S. contributed to conceptualization; S.Y.Z.S. and H.F.S. contributed to data curation; S.Y.Z.S., C.T.N., S.A., W.K.Y., M.O., and H.F.S. contributed to formal analysis; S.Y.Z.S., M.A., M.O., and H.F.S. contributed to funding acquisition; S.Y.Z.S., C.T.N., and S.A. contributed to investigation; S.Y.Z.S., C.T.N., S.A., W.K.Y., and H.F.S. contributed to methodology; S.Y.Z.S., C.T.N., and S.A. contributed to project administration; S.Y.Z.S., C.T.N., S.A., K.T., and M.O. contributed to resources; S.Y.Z.S., C.T.N., and S.A. contributed to software; M.A., K.T., and H.F.S. contributed to supervision; S.Y.Z.S. and H.F.S. contributed to validation; S.Y.Z.S. and H.F.S. contributed to visualization; S.Y.Z.S. contributed to writing – original draft; S.Y.Z.S. and H.F.S. contributed to writing – review and editing.

All data generated or analysed during this study are included in this article and its online supplementary files. Further enquiries can be directed to the corresponding author.