This cross-sectional study was carried out at the First Affiliated Hospital of Zhengzhou University from June 2020 to February 2021. It was approved by the Ethics Committee of the First Affiliated Hospital of Zhengzhou University (This trial is registered on the website ¹ and its unique identifier is approved protocol number: 2021-KY-0176-002) and was conducted according to the Declaration of Helsinki. Fifty-one patients undergoing pars plana vitrectomy were included in this study. All patients were informed of the use of their samples and signed informed consent. In the screening cohort (15 eyes of 15 patients), five patients diagnosed with idiopathic macular holes were included as control group, five PDR patients who had no history of ophthalmic therapy were included as no-treatment group, and five PDR patients who were at the seventh day after intravitreal aflibercept injection were listed as treatment group. In the validation cohort (36 eyes of 36 patients), there were respectively 12 patients in the control group, no-treatment group, and treatment group.

PDR patients were diagnosed according to the 2017 Diabetic Retinopathy Preferred Practice Pattern Guideline issued on the website ² by the American Academy of Ophthalmology: patients with diabetic retinopathy should be considered to be at the PDR stage if they have one of the following manifestations: neovascularization; vitreous or preretinal hemorrhage. In our research, to minimize the effect of blood mixture, only patients with severe fibrous proliferation or tractional retinal detachment were included as PDR group. IMH patients were diagnosed based on clinical manifestation, and fundus and optical coherence tomography results. For the exclusion criteria, patients with dense vitreous hemorrhage, acute or chronic infection, history of ocular treatment, ocular trauma, other ocular diseases, severe heart, liver, and kidney diseases, and infection were excluded.

Under local anesthesia, a 3-port 23-gauge pars plana vitrectomy was performed on all patients by the same surgeon. At the onset of surgery, 1 ml vitreous humor was extracted from the core of the vitreous cavity into a syringe connected to the cutter after the infusion was closed. Before that, the testing water in the cutter tube was pushed out by the syringe. Then the sample was injected into an enzyme-free Eppendorf (EP) tube, put on ice immediately, and centrifuged for 10 min at 12,000×g to remove any cell debris and blood admixture. All samples were stored at −80°C for further use. To minimize bias, all samples were labeled by subject ID according to the randomization protocol, which was concealed in identical sealed envelopes, and the investigators were unaware of the group labels.

Frozen vitreous samples were thawed on ice for an hour prior to centrifugation at 12,000×g for 5 min at 4°C. Total RNA was extracted from 500 μl vitreous using Trizol reagent (Invitrogen, Carlsbad, CA, USA). The purity and concentration of total RNA were quantified by Nano Drop ND-1000 Spectrophotometer (Thermo Fisher Scientific, MA, USA), and an Agilent 2100 bioanalyzer (Thermo Fisher Scientific) was used to assess its integrity.

Total RNA (1 µg) per sample was used to conduct miRNA sequencing library. After size fractionating by 15% denaturing polyacrylamide gel electrophoresis, the RNA mixture was eventually isolated into small RNA fragments with the length of 18–30 nt, the 3′ and 5′ ends of which were ligated with adenylated adapters annealed to unique molecular identifiers (UMIs). Then, the tagged small RNAs were retro-transcribed into single-stranded complementary DNA (cDNA) with Superscript II Reverse Transcriptase (Invitrogen, USA) and cDNAs were amplified by PCR with PCR Primer Cocktail and PCR Mix. The PCR constructs were then separated by agarose gel electrophoresis and target fragments with the size of 110–130 bp were excised and purified by QIAquick Gel Extraction Kit (QIAGEN, Valencia, CA). The final library was qualified and quantified in two methods: the quality and production were tested by real-time quantitative PCR (TaqMan Probe) and the distribution of the fragment size was checked by Agilent 2100 bioanalyzer (Thermo Fisher Scientific). Library preparations were pooled and sequenced by the BGISEQ-500 platform (BGI-Shenzhen, China).

Bioinformatics tools were applied to extract meaningful information from raw sequencing data. After filtering adapter sequences and low-quality tags, clean tags were aligned to the reference genome and miRbase with Bowtie2² (Langmead et al., 2009). Subsequently, Rfam mapping was performed to identify tags that matched rRNA, tRNA, snRNA, snoRNA, and other ncRNAs in cmsearch-tool ³ (Nawrocki and Eddy 2013). The software miRDeep2 ⁴ was used to predict novel miRNAs by exploring the secondary structure (Friedlander et al., 2008). Theoretically, the UMI counts were log-normalized to calculate the absolute expression level of miRNA. Then, the target genes of miRNAs were commonly predicted by the online tools, including RNAhybrid ⁵ (Kruger and Rehmsmeier 2006), miRanda ⁶ (John, Enright et al., 2004), and TargetScan ⁷ (Agarwal, Bell et al., 2015). Differentially expressed miRNAs (DEMIs) were analyzed by DEGseq, ⁸ and only the genes (miRNAs) with the absolute value of Log2Ratio ≥1 and Q value ≤0.01 were considered as DEMIs (Wang, Feng et al., 2010). Gene Ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment were performed to annotate the function of target genes using phyper, a function of R, version 3.1.2. After being corrected by Bonferroni method, p-value ≤0.05 was used as the threshold.

The microarray expression profiling data of total RNA (GSE160310) from 80 human post-mortem retinal samples based on GPL20301 Illumina HiSeq 4000 (Homo sapiens) were downloaded from the Gene Expression Omnibus (GEO) database, ⁹ which included 20, 20, 35, and 5 samples of the four patient groups (healthy control, diabetic, NPDR, PDR), respectively. All data were acquired and applied strictly according to the publication guidelines and data access policies of GEO database. Differential expression genes were analyzed by GEO2R online tool and corrected p-value <0.05 and | logFC | >1.5 were regarded as the threshold.

The validation experiment was conducted with quantitative real-time PCR (qRT-PCR) in another cohort including 12 IMH patients, 12 PDR no-treatment patients, and 12 PDR treatment patients. Total RNA from vitreous sample was extracted according to the manufacturer’s protocol of Trizol reagent (Servicebio, Wuhan, China) with A.C. elegans miR-39, a custom spike-in control. Primer sequences for retrotranscription and qRT-PCR of each miRNA were designed and purchased from Servicebio (Wuhan, China). In this research, cel-miR-39-3p (forward (F), 5′-ACA​CTC​CAG​CTG​GGG​TCA​CCG​GGT​GTA​AAT​C-3′; reverse (R), 5′-CTC​AAC​TGG​TGT​CGT​GGA​GTC​GGC​AAT​TCA​GTT​GAG​CAA​GCT​GA-3′) was used as endogenous control for normalization of the miRNA expression data. Then, the RNA was reverse-transcribed into cDNA using Servicebio RT First Strand cDNA Synthesis Kit (Servicebio, Wuhan, China) and qRT-PCR was performed with 2× SYBR Green qPCR Master Mix (High ROX) (Servicebio) in accordance with the manufacturer’s instructions. Raw qRT-PCR data were analyzed using QuantStudio Real-Time PCR System software (Applied Biosystems). All data represented the average results from three independent experiments. Cycle threshold (Ct) values for each sample were converted to arbitrary absolute amounts (2^−ΔΔCT).

The statistical analysis of data was conducted using IBM SPSS-25 Statistics (IBM SPSS Statistics for Windows, version 20.0; IBM Corp., Armonk, NY, USA). The categorical variables, gender and hypertension, were presented as frequency counts and percent, and were analyzed using Pearson’s χ² or Fisher’s exact test. The continuous variables were summarized as means ± SDs and medians with ranges or quartiles, and one-way ANOVA, Student t-test, or non-parametric testing, including the Kruskal-Wallis test, were conducted for their analysis. A p-value <0.05 indicated statistical significance.

A total of 51 patients were eligible in the study (26 males, 25 females, average age 65.65 ± 5.56), the baseline characteristics of which are presented in Table 1. All raw clinical data was listed in Supplementary File S1, Table 1. As shown in Table 1, there was no significant difference in age or gender distribution between PDR patients and controls in any of the three groups (p > 0.05). In terms of general condition parameters, including hypertension, diabetes duration, and HbA1c, no significant differences of them were observed between all groups in both screening cohort and validation cohort (p > 0.05), except for the comparison of HbA1c between control group and PDR patients.

Total RNAs with good quality were obtained from 15 vitreous samples at an average concentration of 3.70 ± 5.01 ng per 20 µl of sample. In the sequencing data, the average raw tag counts detected from control, no-treatment group, and treatment group were 23,628,974, 19,404,071, and 24,028,428, respectively. After filtrating low-quality tags, the average clean tag counts of miRNA gained from each group were respectively 17,614,269, 16,870,205, and 20,184,021. Also, the average percentage of clean tags for all samples was 82.3%. Under the default local alignment settings, 43.5% clean tags were successfully aligned to reference genome performed by Bowtie2. Ultimately, a total of 405 known miRNAs matched to the miRbase database were identified, and 60 novel miRNAs were predicted by miRDeep2 (Supplementary Material: Table 2). Among the 465 miRNAs, 242, 317, and 334 miRNAs were respectively obtained from each group. A Venn graph was conducted to describe the distribution of miRNAs in each group (Figure 1A).

In our previous research, a detailed analysis about miRNA profiles of control and no-treatment group has been described (Guo, Zhou et al., 2021). To further understand the complex regulation and expression of genes happening in the intraocular environment after intravitreal aflibercept injection, we selected and analyzed the differentially expressed miRNAs among three groups. In the comparison groups of no-treatment/control, treatment/control, and treatment/no-treatment, which were encoded as B/A, C/A, and C/B, respectively, 179 up-DEMIs and 46 down-DEMIs (B/A), 121 up-DEMIs and 51 down-DEMIs (C/A), and 73 up-DEMIs and 131 down-DEMIs (C/B) were screened out (shown in Figure 1B). In total, 29 miRNAs were selected as aflibercept-induced miRNAs, among which 20 miRNAs commonly existed in B/A up-DEMIs, C/A up-DEMIs, and C/B down-DEMIs (shown in Figure 1C); 2 miRNAs commonly existed in B/A down-DEMIs, C/A down-DEMIs, and C/B up-DEMIs (shown in Figure 1D); 5 miRNAs commonly existed in B/A up-DEMIs, C/A down-DEMIs, and C/B down-DEMIs (shown in Figure 1E); and 2 miRNAs commonly existed in B/A down-DEMIs, C/A up-DEMIs, and C/B up-DEMIs (shown in Figure 1F). Moreover, the expression of these 29 aflibercept-induced miRNAs in each group and even per sample is fully described in Figure 2 and Table 2. Samples A1–5, B1–5, and C1–5 in Figure 2 respectively represent five samples in control, no-treatment group, and treatment group.

Based on the prediction results by RNAhybrid, Miranda, and TargetScan, the intersection of 129,936 human genes (including 19,984 mRNAs) were identified as putative targets of these aflibercept-induced miRNAs. To further gain insight into the molecular mechanisms underlying PDR, GO function and KEGG pathways enrichment analysis was performed in 19,984 target mRNAs. On the basis of the GO enrichment analysis, approximately 12,110 biological processes, 1,732 cellular components, and 4,233 molecular functions were identified with Q value (an adjusted p-value) <0.05 as the threshold, and the top 20 terms are shown in Figure 3A–C. The top listed GO terms related to cellular components (CC), molecular functions (MF), and biological process (BP) were respectively membrane (GO:0016020), protein binding (GO:0005515), and phosphorylation (GO:0016310). In the KEGG pathway analysis, approximately 392 pathways were enriched at Q value (adjusted p) <0.05, the top five of which were metabolic pathways (hsa01100; 1,528 mRNA), pathways in cancer (hsa05200; 604 mRNA), biosynthesis of secondary metabolites (hsa01110; 461 mRNA), axon guidance (hsa04360, 196 mRNA), and biosynthesis of antibiotics (hsa01130, 276 mRNA) (Figure 3D). What is more, 67 target mRNAs play roles in the VEGF signaling pathway with a Q value of 0.008 according to the enrichment result (hsa04370), and a network was conducted to illustrate the target relationship between the 29 miRNAs and mRNAs in VEGF signaling pathway (Figure 3E), the specially marked parts of which show seven miRNAs (hsa-miR-3074–5p, hsa-miR-125b-5p, hsa-miR-197-3p, hsa-miR-204-3p, hsa-miR-24-3p, hsa-miR-574-3p, hsa-miR-34a-5p) that directly target VEGFA.

After a comprehensive analysis, three miRNAs were randomly selected among the 29 aflibercept-induced miRNAs to confirm the expression profiles. The validation result was listed in Supplementary File S1, Table 3. As expected, the qRT-PCR result was consistent with the NGS result, indicating the reliability of the expression profiles (Figure 4; Table 3). In the validation cohort, hsa-miR-3184-3p, hsa-miR-24-3p, and hsa-miR-197-3p were significantly upregulated in the no-treatment group compared with the control group, and significantly downregulated in the treatment group compared with the no-treatment group (p < 0.05). There was no significant difference in the expression of has-miR-24-3p and hsa-miR-3184-3p between control group and treatment group (p = 0.078, p = 0.543). Hsa-miR-197-3p was upregulated in the treatment group compared with the control group.

To further understand the expression of 19,984 target mRNAs in PDR patients, the raw data (GSE160310) about microarray expression profiling of total RNA in retinal tissues of human post-mortem donors were retrieved from the GEO database, and the data of 20 healthy controls and 5 PDR donors were included in this study. In the new profiling result, 204 protein-coding mRNAs were found to be differentially expressed in the PDR group compared with healthy control, including 179 upregulated and 25 downregulated mRNAs (Figure 5A), and finally we only identified 179 overlapped mRNAs between the 19,984 target mRNAs in our experiment and 204 deregulated protein-coding mRNAs in the GSE160310 profiling data (Figure 5B). The KEGG pathway annotation for 179 overlapped mRNAs is shown in Figure 5C, and it revealed that immunity (complement and coagulation cascades, hsa04610), infections (pertussis, hsa05133; Staphylococcus aureus infection, hsa05105), inflammation (cytokine–cytokine receptor interaction, hsa05322), immune diseases (systemic lupus erythematosus fiber, hsa05322), formation of fibrous tissue (TGF-beta signaling pathway, hsa04350), cell growth, and death-related pathways (p53 signaling pathway, hsa04115) were highly involved.