This study included 117 participants (M 62/ F 55), recruited from the University Clinical Hospital of Gdańsk, a tertiary reference center in Poland, between March 2021 and June 2022. The cohort comprised 81 patients diagnosed with ADPKD and 15 patients with CKD of other etiologies, including IgA nephropathy or other/unknown primary diseases, excluding diabetes mellitus (DM). The control group consisted of 21 individuals, recruited from the general population of the region, who self-reported as healthy. All participants were adults and underwent screening for eligibility based on predefined inclusion and exclusion criteria specific to their respective groups. Diagnosis of ADPKD was confirmed according to clinical criteria, Ravine-Pei [7].

ADPKD and CKD group’s exclusion criteria: lack of informed consent; diagnosis of DM as an underlying cause of CKD or as a comorbidity, due to anticipated alterations in mineral and bone metabolism; current or previous immunosuppressive therapy, including corticosteroids, owing to their known impact on bone phenotype; estimated glomerular filtration rate (eGFR) < 15 ml/min/1.73 m², indicative of advanced CKD-associated mineral and bone disorder; known allergy to local anesthetic agents required for bone material strength index (BMSi) assessment; presence of leg oedema or obesity, which may interfere with probe penetration and compromise the accuracy of BMSi measurements.

Healthy volunteers’ exclusion criteria: lack of informed consent, known allergy to local anesthetic agents required for BMSi assessment; presence of leg oedema or obesity, which may affect probe penetration and the accuracy of BMSi measurements.

Further stratification of ADPKD participants was based on the stages of CKD according to the Kidney Disease: Improving Global Outcomes (KDIGO) guidelines [8] and the stages of Mayo Clinic Imaging Classification (MIC) [9]. The Chronic Kidney Disease Epidemiology Collaboration (CKD-EPI) equation was used to estimate GFR from measured serum creatinine levels [10].

Class 1 ADPKD patients were stratified into five subclasses based on age and height-adjusted total kidney volume (HtTKV) assessed with magnetic resonance imaging (MRI). According to these parameters, Class 1 subgroups named A, B, C, D and E are identified, and the risk for GFR deterioration is assumed to increase progressively from Class 1 A to Class 1E [9].

Fasting venous blood and urine samples were collected. The biological material was centrifuged immediately and stored (plasma and serum) at − 70 °C until analysis.

Serum samples of creatinine, albumin, calcium, phosphate, magnesium, intact parathyroid hormone (PTH), vitamin D, high-sensitivity C-reactive protein (CRP) and ALP were measured by routine automated methods at the certified Central Clinical Laboratory, University Clinical Center, Gdańsk, Poland.

Bone biomarkers were analyzed at the Department of Clinical Chemistry (Swedac accredited no. 1342), Linköping University Hospital, Linköping, Sweden. All samples were assayed with reagents from the same batch. Serum bioactive sclerostin was measured by enzyme-linked immunosorbent assay (ELISA) (Biomedica, Vienna, Austria), intact fibroblast growth factor-23 (FGF23) by ELISA (Kainos Laboratories, Inc., Tokyo, Japan) and C-terminal FGF23 by ELISA (Biomedica, Vienna, Austria). Both serum BALP and TRACP5b were analyzed with chemiluminescent immunoassays assays on the iSYS fully automated instrument (Immunodiagnostic Systems, Ltd., Boldon, UK). Serum intact PINP was measured with the UniQ radioimmunoassay (Aidian Oy, Espoo, Finland).

MR Siemens Magnetom Aera 1,5T (Siemens Healthineers Global) was used to perform MRI. The protocol based on standard guidelines does not require contrast medium injection [10]. Data were analyzed using Syngo.via MR Siemens software.

Areal bone mineral density (aBMD) was measured by dual-energy X-ray absorptiometry (DXA) in four regions: lumbar spine, proximal femur, non-dominant forearm and whole body. Hologic Discovery Wi bone densitometer was used for this purpose (Hologic, Inc., Marlborough, MA, USA). A standard protocol for positioning was used to reduce aBMD variability. Results were reported according to the International Society for Clinical Densitometry recommendations by one clinician [11]. Lumbar spine scans were re-analyzed using TBS iNsight computer software (Medimaps, Merignac, France version 3.0.2) to calculate trabecular bone score (TBS).

The bone mechanical properties were assessed using OsteoProbe® (Active Life Scientific, Santa Barbara, CA, USA) to measure the BMSi. Following local anesthesia, the handheld OsteoProbe® was inserted through the skin at the midshaft of the right tibia, positioned at the mean distance between the distal apex of the patella and the medial malleolus. The device was advanced until it reached the bone surface, which was indented upon activation.

A minimum of 8 to 10 measurements were performed from a single skin insertion site in a clockwise direction. After completing bone measurements, additional measurements were conducted on a polymethylmethacrylate (PMMA) plastic calibration phantom. The BMSi outcome measure was calculated as 100 times the ratio of the indentation distance in the PMMA standard divided by the indentation distance in bone.

Quantitative and categorical variables were presented as medians with interquartile ranges and numbers with percentages, respectively. Patients were categorized into various groups and compared. The Wilcoxon rank sum test and the chi-squared test were used to identify statistically significant differences between two groups for quantitative and categorical variables, respectively. The Kruskal-Wallis test and the chi-squared test were used to identify statistically significant differences between three groups for quantitative and categorical variables, respectively.

Two complementary heat maps were generated to visualise the distribution of mineral and bone parameters across disease stages within the ADPKD cohort only. The first heat map was based on CKD staging reflecting disease progression according to kidney function (eGFR). The second was constructed according to the MIC, which stratifies patients morphologically based on HtTKV. This approach allowed parallel visualisation of the skeletal phenotype along both functional and structural dimensions of ADPKD progression. It enabled us to identify domain-specific patterns (bone mass, bone strength, biomarkers, bone turnover) and compare their distribution between two classifications. Z-scores were computed separately for each variable, based on its own mean and standard deviation within the analysed cohort. Z-scores are not absolute differences, and are intended to give an overview of patterns rather than to indicate statistical significance. Accordingly, the color intensity in the heat maps represents the relative deviation of a given parameter from its mean, not absolute values. As a result, similar colors may correspond to different absolute Z-scores across variables, reflecting independent scaling of each parameter.

We applied a rank-based linear regression method. Compared to conventional linear regression, rank-based regression is inherently more robust, as outlying observations exert a far weaker influence on the resulting model. Estimation of the regression coefficients was carried out by minimizing the convex and piecewise linear (CPL) regression-rank criterion function, as proposed by Bobrowski [12, 13]. Normalized feature importance scores (0–100%, summing to 100%) were calculated to reflect each variable’s influence. Conceptually, this approach maximizes the probability that, when comparing two patients within the ADPKD cohort, the model will correctly identify which one presents with a more advanced disease stage. The analysis was carried out across three distinct groups designated as ranks 1, 2, and 3, respectively: CKD G1-2, CKD 3a, and CKD G3b-4, according to the CGA; and Mayo A-B, C, and DE, according to the MIC classification.

Statistical analyses were performed using R version 4.5.1 and own implementation in Python 3.9 environment for the ranked regression CPL method. Results were considered statistically significant at a p value < 0.05, unless otherwise indicated.