Genotype–Phenotype Associations in an X-Linked Retinoschisis Patient Cohort: The Molecular Dynamic Insight and a Promising SD-OCT Indicator

Abstract

Purpose: This study investigated a three-dimensional indicator in spectral-domain optical coherence tomography (SD-OCT) and established phenotype-genotype correlation in X-linked retinoschisis (XLRS).

Methods: Thirty-seven patients with XLRS underwent comprehensive ophthalmic examinations, including visual acuity (VA), fundus examination, electroretinogram (ERG), and SD-OCT. SD-OCT parameters of central foveal thickness (CFT), cyst cavity volume (CCV), and photoreceptor outer segment length were assessed. CCV was defined as the sum of the areas of cyst cavities in uential B-scans, measured automatically by self-developed software (OCT–CCSEG). Structural changes of the protein associated with missense variants were quantified by molecular dynamics (MD). The correlation between genotype and phenotype was analyzed.

Results: Twenty-seven different RS1 variants were identified, including a novel variant c.336_337insT(p.L113Sfs*8). The average age of onset was 14.76 ± 15.75 years, and the mean VA was 0.84 ± 0.43 logMAR. The mean CCV was 1.69 ± 1.87 mm³, correlating significantly with CFT (R = 0.66; P < 0.01). In the genotype–phenotype analysis of missense variants, CCV significantly correlated with the structural effect on the protein of mutational changes referred to as wild type, including root-mean-square deviation (R = 0.34; P = 0.04), solvent accessible surface area (R = 0.38; P = 0.02), and surface hydrophobic area (R = 0.37; P = 0.03). The amplitude of scotopic 3.0 ERG a-waves and b-waves significantly correlated with the percentage change of the β-strand in the secondary structure (a-wave: R = −0.58, P < 0.01; b-wave: R = −0.53, P < 0.01).

Conclusions: CCV is a promising indicator to quantify the structural disorganization of XLRS retina. The OCT–CCSEG software calculated CCV automatically, potentially facilitating prognosis assessment and development of personalized treatment. Moreover, MD-involved genotype–phenotype analysis suggests an association between protein structural alterations and XLRS severity measured by CCV and ERG.

X-linked retinoschisis (XLRS; OMIM 312700) is an X-linked recessive inherited retinal disorder with an estimated prevalence ranging from 1 in 5000 to 1 in 25,000.^1–4 Patients with XLRS typically present with visual acuity loss, characterized by bilateral foveal schisis arising from the splitting of the inner retinal layers and a negative electroretinogram (ERG) waveform with significant attenuation of the b-wave amplitude.^2,5–8 XLRS is caused by variants in the RS1 gene, which encodes an adhesive protein, retinoschisin (RS1), which is essential for the structural and functional integrity of the retina. 

The phenotype of XLRS exhibits considerable heterogeneity, presenting practical challenges in disease management. Prior studies have found that cystoid changes in retinoschisis can involve various retinal layers, extending from the retinal nerve fiber layer to the nuclear layer, with the extent of the schisis showing significant variability.^9–11 This highlights the necessity to identify appropriate clinical parameters that accurately reflect the effect of the disease on retinal structural changes. For the assessment of structural retina changes and prognosis of XLRS, spectral-domain optical coherence tomography (SD-OCT) is an effective approach to examine the multiple properties of retinal schisis and provide quantifiable indicators of structural changes. However, existing studies have mainly focused on one- or two-dimensional SD-OCT indicators, such as total retinal thickness, central foveal thickness (CFT), and area of the macular schisis cavity (AMS).^12,13 These indicators have limitations with regard to reflecting the actual structural disorganization. Thus, we were motivated to consider three-dimensional (3D) volumetric parameters in SD-OCT, which potentially offer a more precise evaluation of retinal structural characteristics in patients with XLRS. 

In addition to the diverse degree of structural lesions in the XLRS retina, considerable variability is observed in retinal function among XLRS patients, as manifested primarily in the ERG a-wave, b-wave, and the b-wave/a-wave ratio. To explain the broad clinical heterogeneity of XLRS, several studies have sought to establish the genotype–phenotype correlation of the disease. Pioneer studies mainly aimed to correlate clinical features, such as visual acuity, foveal pathology, complications, and ERG b-wave amplitude, with variant types but failed to establish a clear correlation.^8,14–17 Recent studies have provided insight into the molecular mechanisms involved. Research conducted by Sergeev et al.^18,19 correlated molecular models of RS1 obtained from molecular dynamics (MD) simulations with retinal function and suggested an association between the structural alteration of protein and the severity of XLRS as determined by ERG b-wave/a-wave ratios. In addition, Sudha's research investigated the association between the molecular characterization of the equilibrated structure of the mutant RS1 protein and its secretion profile.²⁰ These studies evaluated the impact of a variant at the atomic level, shedding new light on the genotype–phenotype analysis of XLRS. 

In this study, we aimed to analyze the clinical and genetic characteristics of 37 Chinese XLRS patients. Our findings revealed a promising method for obtaining 3D volumetric parameters using SD-OCT that can be applied along with other indicators to further investigate the correlation between phenotype and genotype in XLRS. Moreover, we gained a novel perspective on the genotype–phenotype correlation analysis of XLRS from molecular dynamic insights. 

Methods

Recruitment of Subjects

Thirty-seven patients with molecularly confirmed XLRS were recruited from the Department of Ophthalmology, Peking Union Medical College Hospital (PUMCH), Beijing, China, and Beijing Mei'ermu Hospital, Beijing, China. Informed consent was collected from them or their legal guardians before enrollment. The PUMCH review board approved this study (no. JS-2059), which was conducted in accordance with the Guidance on Sample Collection of Human Genetic Diseases by the Ministry of Public Health of China, as well as the tenets of the Declaration of Helsinki. 

Clinical Evaluations

Detailed medical and family histories were recorded for each patient. Patients underwent comprehensive ophthalmological examinations, including visual acuity, slit-lamp biomicroscopy, dilated indirect ophthalmoscopy, fundus photography, fundus autofluorescence photography (SPECTRALIS; Heidelberg Engineering, Heidelberg, Germany), SD-OCT (SPECTRALIS; Heidelberg Engineering), and full-field standard ERG (RETIport; an-vision, Salt Lake City, UT, USA). VA was assessed using the Snellen chart and converted to logMAR units. For OCT image acquisition, a skilled technician performed an SD-OCT examination on all of the patients with XLRS. The scanning area was 5.6 × 5.6 mm in in the XY plane, centered at the fovea, generating a set of 25 horizontal cross-sectional B-scan images. Each OCT scan was 5.6 × 1.9 mm in the XZ plane (Fig. 1). ERG was performed with corneal ERG-Jet contact lens electrodes (Fabrinal, Neuchatel, Switzerland) or DTL electrodes (Roland Consult Stasche & Finger GmbH, Brandenburg ad Havel, Germany). The test protocol followed the standards of the International Society for Clinical Electrophysiology of Vision. 

Figure 1.

 

Illustration of the SD-OCT scanning pattern. The scanning pattern area was 5.6 × 5.6 mm in in the XY plane, centered at the fovea, which generated a set of 25 horizontal cross-sectional B-scan images. Each OCT scan was 5.6 × 1.9 mm in the XZ plane.

To further assess SD-OCT imaging biomarkers, the images were processed through ImageJ (National Institutes of Health, Bethesda, MD, USA) to calculate the CFT and photoreceptor outer segment (PROS) length. The PROS length was determined as the perpendicular distance between the posterior surface of the inner segment (IS)/outer segment (OS) junction and the anterior surface of the retinal pigment epithelium (RPE), taken at the foveal center on a horizontal central SD-OCT image (Supplementary Fig. S1). The retinal layers with foveoschisis were documented. To accurately quantify the macular cyst cavities in SD-OCT images, we defined and calculated cyst cavity volume (CCV) as the sum of the areas of cyst cavities across 25 sequential horizontal cross-sectional B-scans, multiplied by the distance between these scans. Mathematically, the CCV can be approximated by:  

\begin{eqnarray*}{\rm{CCV}} = \mathop \int \nolimits_0^L {S_{cyst}}\ dl \approx \mathop \sum \limits_{i = 1}^{25} {S_{cyst\ in\ i - scan}}{\rm{\Delta }}l = {\rm{\Delta }}l\mathop \sum \limits_{i = 1}^{25} {S_{cyst\ in\ i - scan}}\ \ \end{eqnarray*}

where Δl represents the interscan distance. It was measured automatically by a self-developed software system, OCT cyst cavity segmentation (OCT–CCSEG), using Python. Specifically, the system consists of five steps: (1) detecting the retinal layers, inner limiting membrane (ILM), and RPE by OCT segmentation with MATLAB software OCTSEG (MathWorks, Natick, MA, USA)²¹; (2) initially segmenting candidate cyst cavities by identification of the contours of all dark regions between the ILM and RPE using the Chan–Vese algorithm; (3) selecting cyst cavities from the candidates based on their intensity, shape, and position; (4) calculating the cyst cavity area in each scan; and (5) computing CCV from 25 SD-OCT scans. The entire framework is depicted in Figure 2. 

Figure 2.

 

The CCV automated calculation system OCT–CCSEG. The input is a group of OCT scans, and the output is CCV. The system includes five modules: image boundaries detection module, image segmentation module, filter module, cavity area calculation module, and cavity volume calculation module.

Genetic Analysis of RS1 Variants

The procedure and variant filtering were conducted as previously described.²² Targeted panel sequencing of 256 known retinal disease genes and whole exome sequencing were used to identify pathogenic variants. Briefly, the genomic DNA of all probands and their available relatives was isolated from peripheral leukocytes and captured via the SeqCap EZ Choice XL Library (Roche NimbleGen, Pleasanton, CA, USA), then sequenced on the Illumina HiSeq platform (Illumina, San Diego, CA, USA). The raw data were analyzed using NextGENe 2.3.4 software, and the reads were compared to the reference sequence of hg19 from the UCSC Genome Browser (University of California Santa Cruz Genomics Institute, Santa Cruz, CA, USA). A comparison was conducted using the 1000 Genomes Project database, dbSNP database, and Exome Aggregation Consortium database to exclude nonpathogenic polymorphisms. 

Pathogenicity of the missense variants and nucleotide conservation were accessed through the VarCards database.²³ In silico algorithms, including SIFT,²⁴ PolyPhen2,²⁵ MutationTaster,²⁶ MutationAssessor,²⁷ FATHMM,²⁸ PROVEAN,²⁹ VEST3,³⁰ MetaSVM,³¹ MetaLR,³¹ M-CAP,³² CADD,³³ DANN,³⁴ fathmm-MKL,³⁵ GenoCanyon,³⁶ and REVEL³⁷ were applied to predict the pathogenicity of the missense variants. Nucleotide conservation was analyzed using GERP++,³⁸ phastCons,³⁹ phyloP,⁴⁰ and SiPhy.⁴¹ The pathogenicity and conservation of small deletions, insertions, and nonsense variants were assessed by Phylop, PhastCons, and MutationTaster. The potential effects of splicing variants were evaluated by the Human Splicing Finder⁴² and SpliceAI To identify regions in the RS1 protein that are sensitive to mutations, in silico saturation mutagenesis analysis was performed on the wild-type RS1 protein using Rhapsody.⁴³ Putative pathogenic RS1 variants (NM_000330.3) were verified by Sanger sequencing and segregation analysis when family members were available. Variants were classified according to the standards and guidelines of the American College of Medical Genetics and Genomics (ACMG), and the nomenclature of all variants was adjusted to the Human Genome Variation Society guidelines.^44,45 The protein model of human RS1 (PDB-3JD6) was applied to visualize the distribution of the variants within RS1 on the protein structure. The rendering of the protein structure and annotation of interesting regions were achieved by the PyMOL Molecular Graphics System 2.0 (Schrödinger, LLC, New York, NY, USA). 

Molecular Dynamics Simulation and Trajectory Analysis

The structural stability and properties of the wild-type and RS1 mutant proteins were analyzed through MD simulations. We applied the monomeric form of PDB-3JD6⁴⁶ as the initial model for wild-type RS1. The structures of the mutant RS1 were generated by Discovery Studio software (Accelrys, San Diego, CA, USA) with the RS1 wild-type structure as a template. The predicted structures were confirmed through a Ramachandran plot that showed no residues in the disallowed region (Supplementary Fig. S2, Supplementary Table S1). Prior to the MD simulation, we performed an energy minimization for the structures to remove any steric conflicts using the Smart Minimizer tool in Discovery Studio. The MD simulations were executed in four steps, including simulation setup, energy minimization, equilibration, and production simulation, using CHARMM force field and Gromacs 2022.1 software.⁴⁷ The parameters used in the MD simulations are detailed in Supplementary File S1. 

The trajectory files of the MD simulations were processed by GROMACS utilities, including gmx rmsd, gmx sasa, and gmx hbond, to analyze root-mean-square deviation (RMSD), solvent-accessible surface area (SASA), and hydrogen bond. The secondary structure was determined through DSSP⁴⁸ and virtualized in UCSF Chimera.⁴⁹ The hydrophilic and hydrophobic areas were estimated by the software BioLuminate (Schrödinger), and the free energy changes in the protein due to missense variants were evaluated through the FoldX force field.⁵⁰ Alteration in the structural property of a mutant protein in comparison to the wild type was computed and is represented by Δ (delta). For example, the ΔRMSD of a given variant signifies the deviation in RMSD between the mutant protein and the wild type. The figures were generated using the Python Matplotlib 3.7.0 package. 

Statistical Analysis

The statistical analysis in the study was performed by Python. The normality of the data was assessed by the D'Agostino or Pearson algorithm. The t-test was performed to determine the difference in phenotype features between patients and normal controls. Results with P < 0.05 were considered statistically significant. The missing observations were ignored, and analysis was done on the variables present. The time-to-event analysis was performed to analyze the impact of age on VA. The non-numerical VA was transformed into logMAR values, such as hand movements to 2.7 logMAR. Based on the visual impairment criteria of the World Health Organization, we divided VA into four levels: blindness (level 1, VA < 1.30 logMAR), severe visual impairment (level 2, 0.30 logMAR ≤ VA< 1.00 logMAR), moderate visual impairment (level 3, 1.00 logMAR ≤ VA < 0.48 logMAR), and mild visual impairment (level 4, 0.48 logMAR ≤ VA < 0.30 logMAR). The Kaplan–Meier survival analysis was applied to estimate the age at these visual impairment levels, and a spline regression model was used to account for censored data. The ordinal logistic regression was used to explore the annual decline/growth ratio for VA, a-wave amplitude, b-wave amplitude, and the b-wave/a-wave ratio. In the genotype–phenotype analysis, patients were categorized into two groups based on the type of variants they harbored¹⁸: missense variants and null variants (nonsense, insertion, small deletion, splicing, and gross deletions). 

Results

Clinical Characteristics

The clinical characteristics of the 37 Chinese patients are summarized in Table 1, including age of onset, VA, SD-OCT indicators (CFT, CCV, and PROS), and ERG parameters. Detailed clinical information for each enrolled patient is presented in Supplementary Table S2. Representative fundus and SD-OCT images and ERG recordings of patients with XLRS are shown in Figure 3. The mean ± SD age of onset of the 37 patients in our cohort was 14.76 ± 15.75 years (interquartile range [IQR], 5.0–16.0; range, 1.0–58.0), and the mean VA was 0.84 ± 0.43 logMAR (IQR, 0.52–1.3; range, 0.10–2.70). To quantify the structural disorganization of the XLRS retina, SD-OCT was conducted for 34 patients (64 eyes). The mean CFT, PROS, and CCV were calculated to be 0.456 ± 0.235 mm (IQR, 0.317–0.551; range, 0.1–1.321), 21.62 ± 7.13 µm (IQR, 17.69–25.32; range, 0.33–42.58), and 1.69 ± 1.87 mm³ (IQR, 3.36–7.99; range, 0–52.18), respectively. 

Table 1.

 

Clinical Observations of Patients (N = 37)

 

Figure 3.

 

Representative clinical images of XLRS patients. (a) Fundus displaying a spoke-wheel–like schisis at the macula. (b) OCT showing splitting of the inner retinal layers. (c) ERG exhibiting reduced waveforms of rod and cone responses and a negative ERG waveform on scotopic 3.0 ERG.

ERGs were obtained for all patients in this cohort, and the parameters for the scotopic and photopic 3.0 stimuli are shown in Table 1. The scotopic ERG a-wave amplitude, b-wave amplitude, and b-wave/a-wave ratio were significantly lower in our XLRS cohort compared to a control group of 40 normal individuals (P < 0.001) (Fig. 4). It was worth noting that significant correlations were identified between the VA and the a-wave and b-wave amplitudes, as well as the b-wave/a-wave ratio for scotopic ERGs (a-wave: R = −0.30, P = 0.01; b-wave: R = −0.43, P < 0.01; b-wave/a-wave ratio: R = −0.36, P < 0.01). To delve deeper into these associations, we defined the better eye as the eye with better VA. In the results from the better eye, we discovered further significant correlations. Specifically, the better eye VA was significantly correlated with the amplitude of the b-wave (R = −0.39; P = 0.03) and the b-wave/a-wave ratio (R = −0.40; P = 0.03) of the better eye. When gene variants were categorized into two classes (null variants and missense variants), no significant differences were observed between these two groups in terms of ERG parameters and VA. 

Figure 4.

 

Comparison of ERG in patients and the control group. The violin and box plots illustrate scotopic 3.0 a-wave amplitude, scotopic 3.0 b-wave amplitude, and 3.0 scotopic b-wave/a-wave ratio in different groups.

Age Impact on VA

The VA of the patients in the study demonstrated a decreasing trend with increasing age, with an acceleration of this trend observed at older ages (Fig. 5a). Specifically, the trend was characterized by an average annual rate of 0.002 logMAR/y under the age of 10, 0.003 logMAR/y between 10 and 20 years of age, and 0.010 logMAR/y above 20 years old. In the survival analysis, the median ages to reach mild, moderate, and severe impairments were 8 years, 33 years, and 59 years, respectively (Figs. 5b–5d). In the analysis of blindness level, we found that 19.9% of patients were predicted to reach blindness at the age of 60 years (Fig. 5e). 

Figure 5.

 

Age impact on VA. (a) The spline regression model of VA and age. (b–e) Kaplan–Meier survival analysis of different visual function indicators for age: survival curves for best-corrected visual acuity (BCVA) of 0.30 logMAR (b); survival curves for BCVA of 0.48 logMAR (c); survival curves for BCVA of 1.00 logMAR (d); and survival curves for BCVA of 1.30 logMAR (e).

Properties and Qualification of Cysts

The retinal layers with foveoschisis were summarized from the SD-OCT images of 64 eyes. Cysts were observed in the ganglion cell layer (GCL), inner plexiform layer (IPL), inner nuclear layer (INL), outer plexiform layer (OPL), and outer nuclear layer (ONL) with prevalences of 10.9%, 1.6%, 100%, 21.9%, and 31.2%, respectively. Pearson's correlation analysis was conducted to investigate the association among the OCT indicators CFT, CCV, and PROS. In the entire cohort, we found that CCV significantly correlated with CFT (R = 0.307; P = 0.014) and PROS (R = −0.25; P = 0.045). Our data indicated that CCV was not correlated with the b-wave/a-wave ratio in scotopic ERGs (R = −0.016; P = 0.900). Additionally, no statistical difference was found in CCV, CFT, and PROS between the patients with missense variants and those with null variants. The distribution of CCV versus age is shown in Supplementary Figure S3. 

Genetic Characteristics

In this study, 37 patients with genetically confirmed XLRS were enrolled. A total of 27 different pathogenic variants were identified, including one novel variant, c.336_337insT(p.L113Sfs*8). The variants were categorized as 12 missense, four nonsense, four gross deletions, three small deletions, three splicings, and one insertion variants. Most of the variants were located in the discoidin domain (21/27), followed by the RS1 domain (4/27) and the N-terminal leader sequence (2/27). Table 2 presents a comprehensive list of the variants, their associated protein changes, variant types, and protein domains. The genomic location of the variants in the protein structure model of RS1 (PDB-3JD6) and the distribution of the variants among protein domain structure are illustrated in Figure 6. The most common variants were c.637C>T (four patients, 10.8%) and c.214G>A (four patients, 10.8%), followed by c.336_337insT (two patients, 5.2%), c.451T>C (two patients, 5.2%), c.590G>A (two patients, 5.2%), and c.608G>A (two patients, 5.2%). The remaining 21 variants were identified in one patient each. Most variants were located in exons 4, 5, and 6. 

Table 2.

 

Genotype Summary of Patients

 

Figure 6.

 

Distribution of reported variants within RS1. (a) The genomic location of the variants in the homology model of RS1 (PDB-3JD6). The RS1 domain is shown in cyan, the DS domain in lemon, and the C-terminal in pink. (b) The spectrum of the 27 RS1 variants. (c) The distribution of the variants; the frequency of variants in each exon is indicated in the bar plot. The locations of variants in the gene and the protein are marked by lines.

In Silico Analysis

In silico analysis was performed on missense variants, small deletion variants, insertion variants, nonsense variants, and splicing variants in our cohort using various algorithms. The majority of the in silico algorithms, such as PolyPhen2, FATHMM, PROVEAN, VEST3, MetaSVM, MetaLR, M-CAP, CADD, DANN, fathmm-MKL, GenoCanyon, and REVEL, indicated that the 12 missense variants in our cohort were damaging or probably damaging to protein structure and could cause disease. The evaluation of the conservation of these missense variants was performed using GERP, phyloP, PhastCons, and SiPhy, which confirmed their conservation. Evaluation of small deletions, insertions (including the novel variant c.336_337insT), and nonsense variants was performed using phyloP, phastCons, and MutationTaster, revealing that these variants were conserved by phyloP and phastCons and predicted to be disease causing by MutationTaster. The HSF Pro platform (HSF and MaxEnt) and SpliceAI were employed to evaluate the potential effects of the splicing variants. MaxEnt and HSF analysis indicated that variant c.53-1G>A may cause a broken wild-type acceptor site. SpliceAI predicted that this variant could cause acceptor loss within the pre-mRNA position of −1 bp, with a probability of 0.69. Similarly, variant c.327-2A>G was predicted to cause a broken wild-type acceptor site by MaxEnt and HSF. SpliceAI analysis suggested that this variant may cause acceptor loss within the pre-mRNA position of −2 bp, with a high probability of 0.99. For additional information on the results of the in silico analysis, refer to Supplementary Tables S3 to S5. Remarkably, Vijayasarathy et al.⁵¹ previously discussed the variant c.52G>A, which involved the terminal nucleotide of exon 1 and was within the 5′ donor splice site of intron 1. The c.52G>A variant disrupted base pairing between the splice site and small nuclear RNA, thereby reducing normal splicing efficiency. 

To investigate the hot region of variants in RS1, we applied the in silico algorithm Rhapsody, which has the benefit of predicting the effect of substituting any amino acid residue in a protein with any of the other 19 amino acids, enabling in silico site-directed mutagenesis analysis. Consequently, we identified residue changes that may be poorly tolerated, which can provide valuable insights for future research. Through protein-wide, site-directed mutagenesis, we identified the top-scoring variants in RS1 with a Rhapsody score of at least 0.86. Our analysis revealed that the highest scoring variants were located in eight locations in the protein, indicating that these crucial residues have the potential to cause disease when mutated (Supplementary Fig. S4). 

Molecular Dynamics Analysis

Trajectory data obtained from MD simulations were utilized to analyze the impact of missense variants on the RS1 protein structure. The stability of the protein was evaluated by the RMSDs of the backbone atoms, which determined the extent of their alteration or displacement during the simulation. We computed RMSDs for both the wild-type and mutant RS1. The RMSD of the wild-type protein was 0.215 Å. The variants E72K (0.260 Å), R102W (0.218 Å), S134P (0.248 Å), Y151H (0.233 Å), G171E (0.218 Å), R182C (0.232 Å), R197H (0.231 Å), R200C (0.223 Å), P203L (0.218 Å), R213W (0.226 Å), and R213Q (0.230 Å) showed higher deviations compared to the wild type, suggesting a decrease in stability due to increased backbone flexibility. The variant Q80H (0.208 Å) represented the lower deviations, indicating increased rigidity in the backbone. These RMSDs of mutants measure the effect of variants on the structure of the protein. The 3D structures of the wild-type and mutant proteins of RS1 in the MD simulation are displayed in Figure 7, and the RMSD and ΔRMSD values for both wild type and variants are provided in Table 2. 

Figure 7.

 

Structure and RMSD graph of wild-type and variant forms of RS1 in MD simulation. In the 3D structure, helices are shown in red, strand in yellow, coli in cyan, disulfide bonds in orange, and variant spots in magenta. The RMSD graph shows the deviation of the backbone in the wild-type (blue) and variant (red) forms of RS1 during the MD simulation. Their difference is represented in the red box.

An analysis was performed to assess the impact of the variants on the secondary structure of the protein during the MD simulation. An increase in the proportion of α-helical structures was found in all missense variants compared to the wild type. Except for R213Q, the proportion of β-strands was decreased in all missense variants. The percentages of α-helical strands and β-strands in the secondary structure, as well as their changes compared with the wild type, are listed in Table 3. Notably, all missense variants maintained two intact disulfide bonds, which are crucial for protein folding and stability (Table 3). These changes in secondary structure impact the overall stability and function of the protein. 

Table 3.

 

Structure Properties of Variant Proteins

 

Additionally, the chemical properties of the mutant protein were assessed to gain further insight into their impact. The parameters that were calculated included the number of hydrogen bonds within the protein and between the protein and solvent, SASA, and hydrophilic and hydrophobic areas of the protein and their differences compared to the wild type (Table 3). We found that the variants had a greater impact on the hydrophobic area than the hydrophilic area. Changes in the hydrophobic and hydrophilic areas of the protein surface can affect its structure and function, including cell adhesion, secretion, and other intermolecular interactions. This could potentially result in more severe clinical symptoms for patients with vision problems. 

Genotype and Phenotype Correlation for Missense Variants

The correlation matrices and P values for the phenotype and genotype for patients with missense variants are shown in Figure 8. The analysis revealed that the CCV exhibited significant correlations with the structural effect on the protein of mutational changes referred to as wild type, including ΔRMSD (R = 0.34; P = 0.04), ΔSASA (R = 0.38; P = 0.02), and Δ surface hydrophobic area (R = 0.37; P = 0.03). Conversely, CFT showed weaker correlations or non-significant relationships with these three genotype features with correlation coefficients of 0.340 (P = 0.046), 0.383 (P = 0.023), and 0.371 (P = 0.028), respectively. Also, the amplitudes of scotopic 3.0 ERG a-waves and b-waves were significantly correlated with the percentage change of the β strand in the secondary structure (a-wave Δβ-strand: R = −0.58, P < 0.01; b-wave Δβ-strand: R = −0.53, P < 0.01). However, the scotopic 3.0 ERG a-wave/b-wave ratio was not correlated with the Δβ-strand (R = 0.03; P = 0.85). 

Figure 8.

 

Correlation matrix of genotype and phenotype. The value in each cell is the Pearson's correlation coefficient of two variables in the corresponding row and column. The cells with purple circles are significant results (P < 0.05).

Discussion

This study identified the clinical and genetic features of an XLRS patient cohort and established a link between the phenotype and genotype in XLRS from the molecular dynamic insight. The phenotype of XLRS patients is known for its heterogeneity, which was also significant in our cohort. In terms of the structural disorganization observed from SD-OCT images, cysts were found in the GCL, IPL, INL, OPL, and ONL with prevalences of 10.9%, 1.6%, 100%, 21.9%, and 31.2%, respectively. Moreover, the degree of schisis differed extensively among different patients. This heterogeneity posed a challenge for disease management; thus, searching for the appropriate clinical parameter to reflect the severity of effect of the disease on structural changes of the retina is worth investigating. In this study, we explored a potential OCT indicator to quantify the structural disorganization of XLRS and further applied it to compare the retinal structural changes in patients with XLRS with different genotypes. 

Previous OCT studies of XLRS focused on CFT, but this indicator is limited to the depiction of a single OCT scan at a one-dimensional level.^12,13 Lin et al.⁵² proposed utilizing the AMS as a parameter in OCT and explored its correlation with VA in patients with XLRS. The measurement of the AMS was performed through manual annotation and selection, a process that is time consuming and labor intensive. Also, evaluating the AMS remained at the two-dimensional level when describing the cyst cavity and was limited as far as reproducing the real cyst cavity. To quantify the 3D structure of the cyst cavity directly, we defined the total CCV as the sum of the total area of the 25 sequential B-scans of the OCT. However, a challenge in its calculation is the heavy workload of manually annotating images. To the best of our knowledge, there is still no commercial and open-source software packages that support the CCV automated computing. Given the current software vacancies, we self-developed an automated CCV computing platform OCT–CCSEG. 

The OCT–CCSEG platform involves three major steps. Initially, we used the open-source program OCTSEG²¹ to detect the boundaries of ILM and RPE in OCT. The feasibility and open-source nature of this program made it a suitable choice. Next, we applied the Chan–Vese algorithm, which is commonly used in medical imaging to segment the brain, heart, and trachea,^53–55 to identify the dark region between the two boundaries. The robustness of this algorithm against noise was particularly beneficial for detecting cysts. Last, we employed a filter based on regional features such as intensity, shape, and position to eliminate non-cyst regions, which further improves the accuracy of our system. We have validated the accuracy, flexibility, and stability of our system in OCT images for human patients and mouse models (data not shown). The results of our experiments reinforced the robustness of our approach and instilled confidence in its reliable application. Our system also automates all processes, reducing manual labor and increasing efficiency in data analysis. Moreover, the results obtained from the 3D measurements are more accurate and provide a more precise depiction of changes in the cyst cavity compared to other existing OCT biomarkers. 

Furthermore, the OCT–CCSEG system demonstrates versatility in its potential functions. Current natural history and gene therapy studies of XLRS lack clinical indicators to understand and quantify the disease progression and prognosis. Previous OCT studies of XLRS focused on total retinal thickness. However, retinal thickness can be influenced by two factors: thinning caused by the outer retinal atrophy or thickening caused by the presence of cystic cavities. When the total retinal thickness is reduced, it is difficult to determine whether this is due to a reduction in cyst cavity volume or photoreceptor loss. To monitor the disease progression in structural changes, CCV could serve as a more precise indicator. Moreover, the automation and rapid calculation capabilities of our platform facilitate frequent and efficient assessments of changes in CCV throughout disease progression. Considering the significant association between CCV and clinical features of genotype in our patients, our platform may assist in selecting therapeutic windows and endpoints in XLRS intervention trials. 

In addition to its application in XLRS, OCT–CCSEG has potential applications in other retinal diseases due to the high demand for retinal volume indicators. This platform can measure retinal volume in OCT. Specifically, the initial step of the platform entails detecting the ILM and RPE layers in OCT. By calculating the region between these two boundaries, the total volume can be determined, and the difference between this total volume and the cyst volume can accurately measure the volume of the retina. The precise measurements of retinal volume render it an invaluable diagnostic tool for diseases that affect retinal volume. Furthermore, this platform holds potential for application in diseases that impact the microstructure of the retina, such as exudates, cysts, and fluid. These diseases include age-related macular degeneration, retinal vein occlusion, diabetic macular edema, ocular inflammation, and diabetic retinopathy. 

In addition to clinical heterogeneity, gene heterogeneity is highly evident in our cohort. A total of 27 different pathogenic variants were confirmed, including one novel variant. These variants were classified into 12 missense, four nonsense, four gross deletions, three small deletions, three splicings, and one insertion variants. Existing literature mostly agrees on the absence of a clear genotype–phenotype correlation in XLRS. Therefore, it is impossible to draw prognostic conclusions in newly diagnosed patients based solely on the molecular result. Previous studies compared missense variants and null variants with respect to various indicators, including ERG indicators (a-wave amplitude, b-wave amplitude, and b-wave/a-wave ratio), CFT and PROS in OCT, and VA.^16,56,57 However, our study did not find any significant differences between missense variants and null variants in terms of the ERG indicators CFT and PROS, CCV in OCT, or VA. This suggests that simply categorizing the variants as missense and null may not be sufficient to fully describe their impact on the phenotype. Inspired by studies that approached the problem from the perspective of protein structure change,^18–20 we conducted MD simulations for all missense variants to gain a wealth of structural information about the mutant protein, such as RMSD, SASA, surface hydrophobic area, and percentage of β-strand in the secondary structure. These MD metrics measure the impact of variants on protein structure from different perspectives, providing a comprehensive representation of the genotype. For the phenotype, we also included PROS length in our analysis, which has been found to have a strong correlation with VA in patients with XLRS.^56,57 

Our genotype–phenotype analysis of missense variants revealed that CCV was significantly correlated with mutation-induced structural alterations, such as ΔRMSD, ΔSASA, and Δ hydrophobic area. In contrast, the widely used metric of CFT showed a weaker correlation or failed to show significance with these structural alterations. This finding suggests that CCV could potentially serve as an alternative indicator for OCT in genotype–phenotype correlation analyses. Although PROS length was found to have a strong correlation with VA in earlier studies,^56,57 our analysis found that the correlations of PROS length with most genotype and phenotype features were not significant. Furthermore, the correlation between the amplitude of the scotopic a-wave and b-wave with the percentage change of the β-strand in the secondary structure suggested that a change in the secondary structure in the protein might result in a reduction of the amplitude of a-waves and the b-waves in scotopic ERGs. 

Conclusions

We analyzed the clinical and genetic features of 37 Chinese patients with XLRS. In this process, we proposed that CCV could be a promising indicator to quantify the structural disorganization of the XLRS retina in three dimensions. To automate the calculation of CCV, we designed a self-developed software system, OCT–CCSEG. This system was applied to compare the retinal structural changes of patients with different genotypes in our cohort. The results from the MD-involved genotype–phenotype analysis of missense variants implied an association between the predicted protein structural alteration and the severity of XLRS as measured by the CCV and ERG. Moreover, CCV has the potential to serve as a valuable indicator of disease progression and may provide a foundation for determining appropriate therapeutic windows and endpoints in XLRS intervention trials. 

Acknowledgments

The authors thank all of the patients and their family members for participating in the study. 

Supported by grants from Beijing National Natural Science Foundation (L222083), CAMS Innovation Fund for Medical Sciences (CIFMS 2022-I2M-C&T-B-024), Fundamental Research Funds for the Central Universities (3332022007) and National Natural Science Foundation of China (82201192). 

Disclosure: X. Wei, None; H. Li, None; T. Zhu, None; Z. Sun, None; R. Sui, None