Patient enrollment

This prospective pilot study was conducted at the Cornea and Ocular Surface Unit of the San Raffaele Scientific Institute, Milan, Italy. The study was carried out in accordance with the tenets of the Declaration of Helsinki, and it was approved by the Institutional Review Board/Ethics Committee (Comitato Etico Istituto Scientifico Ospedale San Raffaele). Written informed consent was obtained from the patients at enrollment. A schematic representation of the study design is shown in Figure 1 .

Schematic representation of the study design. Patients diagnosed with KC and presenting Vit D insufficiency (<30 ng/mL) were enrolled. They were followed up for 12 months. Full ophthalmic visits were scheduled at months 0, 2, 4, 6, 9, and 12. Serum Vit D and calcium levels were evaluated at months 0, 4, 6, 12. In addition, peripheral blood samples were collected to analyze other serum biomarkers at months 0 and 6.

Briefly, after the screening visit (month 0 [M0]), patients were followed up for 12 months. Follow-up visits were scheduled at 2, 4, 6, 9, and 12 months from M0. A full ophthalmic exam was performed at each visit, which included measures of best spectacle-corrected visual acuity (BSCVA) and corneal tomography and topography. Maximal keratometry (Kmax), thinnest corneal thickness (TCT), and epithelial and stromal minimal thickness were evaluated. At each visit, patients were reminded not to rub. The examination was always performed at the same time of the day (14-15 hours). Additionally, blood samples were collected at the screening visit, and at months 4, 6, and 12.

Patients presenting KC and Vit D insufficiency, that is serum levels less than 30 ng/mL were included in the study. Exclusion criteria were the following: prior surgical procedures of the cornea, including CXL; diagnosis of end-stage KC defined as corneal thickness less than 300 μm and/or extensive apical leucoma and/or corneal hydrops; diagnosis of active keratitis/conjunctivitis. The diagnosis of KC was confirmed by corneal tomography/topography (MS-39 AS-OCT; CSO). At the screening visit, medical history was collected, including the presence of allergies and rubbing habits. All patients were instructed not to rub regardless of having the habit or not. A full ophthalmic exam was performed, including measurement of (1) BSCVA, (2) uncorrected visual acuity (3) Kmax, and (4) epithelial and stromal minimal thickness and TCT by corneal tomography and topography. Patients were treated as per standard of care in case of Vit D insufficiency with oral cholecalciferol supplementation (50,000 IU once a week for the first 3 months). Maintenance treatment was then continued with 50,000 IU once a month up to month 6 (M6).

When CXL was performed during the study, the treated eye (s) was/were not further considered for the analysis. However, if the eye had worsened more than 1 D before CXL, it was included among worsened eyes; while if the eye had worsened less than 1 D before the procedure, it was excluded from further analysis.

Primary outcome

The primary outcome measure was defined as the percentage of patients with Kmax progression of less than 1 D throughout the 12-month follow-up time.

Secondary outcomes

Secondary outcome measures were changes in BSCVA, Kmax, and TCT, as well as Vit D and calcium serum levels, during the 12-month follow-up time.

Clinical parameters

BSCVA, recorded in Snellen equivalents during visits, was converted into logMAR scale values “Counting finger” was converted to 2.0, and “hand motion” was converted to 3.0 logMAR values. ^, BSCVA rate was defined as the difference between BSCVA at a certain month of follow-up and BSCVA at enrolment, divided by the number of months that intercurred between the two-time points, eg, BSCVA rate at M6 = (M6 logMAR—M0 logMAR)/6.

Kmax, TCT, epithelial thickness, and stromal thickness were measured using the MS-39 AS-OCT corneal tomographer/topographer. Kmax rate and TCT rate were calculated in a similar manner as the BSCVA rate.

Thus, an increase in each rate was associated with KC-related worsening of the index parameter: increased BSCVA rates meant that BSCVA had decreased over time, increased Kmax rates meant that Kmax had increased, and increased TCT rates meant that TCT had reduced.

Bioinformatic analysis

FASTQ sequencing reads were adaptor-trimmed and quality-filtered with Trimmomatic, prior to mapping to the hg38 human reference genome ( https://www.gencodegenes.org/human/ ) with STAR. Gene counts were obtained using featureCounts and Genecode basic annotation v31. Downstream gene expression analysis has been done using the R/Bioconductor framework. Sequencing coverage was considered inadequate for 3 samples, which were removed from the analysis. Only genes with a mean of at least 10 read counts per sample were considered for further analyses. Following this criterion, 12,891 protein-coding and 6675 noncoding genes were identified, and for the sake of this analysis, only protein-coding genes were considered.

Principal component analysis

Principal component analysis (PCA) was performed on the 1000 coding genes with the highest amount of variability (ie, explanatory power) in our dataset. Based on this data, the two main dimensions, dimension 1 (26.2% variance) and 2 (12.4% variance) were subsequently taken into consideration, collectively explaining 38.6% of our dataset variance.

UMAP visualization and unbiased clustering analysis

The top 9 components of the PCA, collectively explaining 78.8% of our dataset variance, have been used to compute the 6 nearest neighbors for each sample and generate the Uniform Manifold Approximation and Projection (UMAP) plot via the umap R package and the unbiased clustering using the Leiden algorithm (resolution 0.8) via the igraph R package.

Differential gene expression analysis

Differential gene expression (DGE) analysis was performed in R (version 4.1.1) using DESeq2. Adjusted P values were obtained based on the Benjamini-Hochberg method. Genes with a false discovery rate <10% and absolute log2FoldChange >0.25 were considered significant. Differentially expressed genes were subjected to Gene Ontology analysis using the R/Bioconductor package clusterProfiler.

Variant calling

Variant calling analysis was performed following GATK guidelines for RNA-seq. In particular, STAR’s two-pass mode alignment was performed to improve alignments around novel splice junctions, duplicate reads were identified, RNA alignments at intronic positions were converted to DNA conventions (SplitNCigarReads), and systematic errors in the base quality scores were detected and corrected with the BaseRecalibrator tool. Finally, HaplotypeCaller was used to call variants toward the reference genome. Variants were then hard filtered using VariantFiltration to exclude likely false-positive ones. Only variants passing filter were kept in the downstream analysis.

Sample size

Published data suggest that KC progression of at least 1 D at the apex occurs in more than 80% of cases. We considered a reduction of KC progression to 50% after 12 months to be clinically significant in the Vit D—supplemented group as opposed to 80% in the untreated population. Hence, considering an alpha error of 0.05 and beta error of 0.2, with 0.8 power, 32 patients were necessary. Considering a potential loss at follow-up of 20%, we enrolled 40 patients.

Analysis

For the primary endpoint, we considered each patient separately. For all the other analyses, we considered each eye from the same patient separately. Data were expressed as standard error of the mean. The statistical significance of the differences between the two groups for continuous variables was assessed using the Wilcoxon signed-rank test. For categorical variables, Fisher’s exact test was employed. A two-sided P <SPAN role=presentation tabIndex=0 id=MathJax-Element-1-Frame class=MathJax style="POSITION: relative" data-mathml='≤’>≤≤
≤
.05 was considered statistically significant. RStudio software (RStudio, PBC) and Prism 10.0 (GraphPad Software) were used for the statistical calculations.