Pathogenesis Study Based on High-Throughput Single-Cell Sequencing Analysis Reveals Novel Transcriptional Landscape and Heterogeneity of Retinal Cells in Type 2 Diabetic Mice

All animal studies were approved by the Shanghai Jiao Tong University School of Medicine Animal Care and Use Committee and performed according to the guidelines of the Association for Research in Vision and Ophthalmology in the Statement for the Use of Animals in Ophthalmic and Visual Research. Male homozygous db/db mice (BKS.Cg-Dock7m^+/+ Leprdb/J) and db/m (Dock7m^+/+ Leprdb) heterozygotes from the same colony were purchased from The Jackson Laboratory (Bar Harbor, ME) and housed in a pathogen-free environment with a 12-h light/dark cycle at 22 ± 1°C with free access to food and water.

After the mice were fasted for 6 h, body weights and tail blood glucose levels were monitored monthly from 2 to 8 months of age. For glucose tolerance tests, mice were fasted overnight before being administered 2 g/kg body weight glucose i.p., and tail blood glucose levels were measured 0, 15, 30, 60, 90, and 120 min later with a glucometer. The area under the curve was calculated at 2 and 8 months of age.

Three db/m and db/db mice each were sacrificed at 8 months of age. One eye from each mouse was used for scRNA-seq, and the other was used for generation of retinal cryosections. The globes were removed and dissected in DMEM (Thermo Fisher Scientific, Waltham, MA) on ice. Retinal tissues were dissociated with use of the Worthington Papain Dissociation System (Worthington, Lakewood, NJ). Briefly, retinal samples were shredded with scalpels and digested in Earle’s balanced salt solution containing papain (20 units/mL) and DNAse I (2,000 units/mL) for 20 min at 37°C. Dissociation was stopped with ovomucoid. Intact cells were separated on a single-step discontinuous density gradient and resuspended in Earle’s balanced salt solution with 10% FBS.

Retinal cell suspensions were loaded on the Chromium Single Cell Controller Instrument (10x Genomics) for generation of single-cell gel beads in emulsion. For scRNA-seq, libraries were prepared with a Chromium Single Cell 5′ Reagent Kit, version 3, according to the manufacturer’s protocol. Sequencing was performed on an Illumina HiSeq X Ten System and 150–base pair paired-end reads were generated. One full lane was used per sample, with approximately >90% sequencing saturation.

The Cell Ranger software pipeline (version 3.0.0) from 10x Genomics was used to demultiplex cellular barcodes, map reads to the genome and transcriptome (with use of the STAR aligner), and down-sample reads as required for generation of normalized aggregate data across samples, to produce a matrix of gene counts versus cells. We processed the unique molecular identifier count matrix using the R package Seurat (version 2.3.4). To remove low-quality cells and likely multiple captures (a major concern in microdroplet-based experiments), we retained the cells with >200 genes and <8,000 genes; >400 unique molecular identifiers and <20% mitochondrial RNA. Library size normalization on the filtered matrix was performed in Seurat to obtain normalized counts.

The top variable genes between single cells were identified with the method described by Macosko et al. (12). Briefly, principal components analysis was performed with the variably expressed genes, and the top 25 principal components were applied for cell clustering with use of a graph-based clustering approach at a resolution of 0.2. We performed t-distributed stochastic neighbor embedding (t-SNE) with the same number of principal components and default parameters to visualize the clustering results. These steps were performed in R with the Seurat package. Comparison with canonic cell type markers indicated that the 17 clusters identified corresponded to 11 cell types.

Single nucleotide polymorphism–trait associations were downloaded on 9 December 2020 from the NHGRI-EBI GWAS Catalog for the traits (13). Heat maps were generated with R. Calculated module scores were calculated for gene expression programs in cell types with use of the AddModuleScore function in Seurat.

To quantify alterations in cell-cell communication in the retina during diabetes, we used CellPhoneDB (v2.0) to identify biologically relevant interacting ligand-receptor partners in the scRNA-seq data (14). We defined a ligand or receptor as “expressed” in a particular cell type if 10% of the cells of that type had nonzero read counts for the gene encoding the ligand/receptor. To define cell-cell communication networks, we linked any two cell types where the ligand and receptor were expressed in the former and latter cell types, respectively. For each receptor-ligand pair in the cell-cell communication, we computed the Spearman correlation coefficient between the mean log₂(normalized data + 1) ligand gene expression in the ligand-expressing cell and the logit-transformed proportions of the receptor-expressing cell across samples.

Differentially expressed genes (DEGs) were identified with use of the FindMarkers function (test.use = bimod) in Seurat. P value <0.05 and |log₂fold change| > 0.26 was set as the threshold for significantly differential expression. Heat maps, volcano plots, and violin plots were generated with R. DEGs were enriched for their involvement in various biological pathways with use of KEGG (Kyoto Encyclopedia of Genes and Genomes) Pathway Enrichment.

Müller glia-specific viruses were constructed with the Müller glia-specific adeno-associated virus (AAV) serotype shH10Y and a ubiquitous synthetic cytomegalovirus promoter (15). AAV vectors expressing either retinaldehyde-binding protein 1 (RLBP1) (ShH10Y.RLBP1) or GFP (ShH10Y.GFP) were used. The animals were randomly divided into four groups: untreated db/m mice, untreated db/db mice, db/db mice treated with intravitreal ShH10Y.RLBP1, and db/db mice treated with intravitreal ShH10Y.GFP. Intravitreal injection of corresponding virus with 4 μL/eye was performed in the mice at 4 months of age. Untreated mice received equal amounts of PBS. The retinas were harvested at 8 months of age. For determination of the transfection efficiency, cryosections were generated from the retinas of db/db mice injected with ShH10Y.GFP.

Immunofluorescence, Western blot, ELISA, electroretinogram, TUNEL assays, trypsin digest for retinal vascular architecture, and detection of vascular permeability are described in Supplementary Materials.

Data were analyzed with GraphPad Prism 8.4.0 (GraphPad Software, San Diego, CA) and are presented as the mean ± SD. Normality was analyzed with the Kolmogorov-Smirnov test. Statistical differences between two groups were compared with unpaired two-tailed Student t test, and those between multiple groups were examined with ANOVA with a Bonferroni correction. P < 0.05 was considered statistically significant.

The data sets have been deposited into National Center for Biotechnology Information Sequence Read Archive (SRA) database under accession number PRJNA653629. All other data generated during the current study are available from the corresponding author upon reasonable request.

The experimental set consisted of retina cells from three control (db/m) and three diabetic (db/db) mice. Their fasting body weights and blood glucose levels were measured monthly from 2 to 8 months of age (Supplementary Table 1). All mice underwent glucose tolerance tests at 2 and 8 months of age, and their glycemic responses at these time points are shown in Supplementary Fig. 1A and B. At 2 months of age, db/db mice already displayed glucose intolerance, as shown by the increased area under the glycemic curve in comparison with that for control mice (P < 0.001 [Supplementary Fig. 1C]). The tendency for glucose intolerance in db/db mice continued to increase until 8 months of age (P < 0.001 [Supplementary Fig. 1C]). Isolated retinas from control and diabetic mice were dissociated into single-cell suspensions. After scRNA-seq and quality control filters, 34,010 and 51,558 retinal cells derived from control and diabetic mice were sequenced, respectively. The cells were classified into 17 transcriptionally distinct clusters (Fig. 1A). Based on the expression of well-established markers, we determined cell identifiers such as rods [rhodopsin (Rho) and G protein subunit α transducin 1 (Gnat1)], cone bipolar cells [secretagogin (Scgn)], rod bipolar cells [Sebox], Müller glia [glutamate-ammonia ligase (Glul) and RLBP1], amacrine cells [solute carrier family 32 member 1 (Slc32a1) and glutamate decarboxylase 1 (Gad1)], retinal ganglion cells [RGCs] [POU class 4 homeobox 1 (Pou4f1) and synuclein γ (Sncg)], cones [opsin 1 (cone pigments) medium-wave-sensitive (color blindness, deutan) (Opn1mw) and opsin 1 (cone pigments) short-wave-sensitive (color blindness, tritan) (Opn1sw)], microglia [allograft inflammatory factor 1 (Aif1) and C-X3-C motif chemokine receptor 1 (Cx3cr1)], vascular endothelial cells [cadherin 5 (Cdh5) and platelet and endothelial cell adhesion molecule 1 (Pecam1)], pericytes [potassium voltage-gated channel subfamily J member 8 (Kcnj8) and platelet-derived growth factor receptor β (Pdgfrb)], and horizontal cells [LIM homeobox 1 (Lhx1)] (16,17) (Fig. 1B and C and Supplementary Fig. 3). The cell clusters were not due to batch or technical effects (Supplementary Fig. 2).

To determine which retinal cell types are vulnerable to DR, we compared the proportions of various cell clusters in the retinas of diabetic and control mice (Fig. 1D). Moreover, genetic variants have been reproducibly linked to DR through population-based GWAS (13). To uncover potential links between the identified cell clusters and DR and PDR, we quantified the overlap between candidate loci from GWAS and the expression signatures of the retinal cell clusters in db/m and db/db mice (Fig. 2A and B). For example, acid-sensing ion channel 5 expression increased in cluster 3 but decreased in clusters 1 and 9. In summary, our single retinal cell sequencing data set can be used to identify the cell types relevant to DR and PDR (Supplementary Fig. 4). Rods, rod bipolar cells, cones, and vascular endothelial cells were the most intimate relative of PDR risk. Furthermore, among Müller glia clusters, cluster 6 was more strongly linked to retinal changes between db/m and db/db mice than the others.

To diagram the wiring of cell-cell interactions more generally between diabetic and control mice, we mapped receptor-ligand pairs onto cell types to construct a putative cell-cell communication across disease states (Fig. 3A). The most active cell types in both states were RGCs, vascular endothelial cells, and Müller glia. Furthermore, microglial activity was increased in the retinal cell-cell communication networks of db/db mice in comparison with those of db/m mice. Conversely, other cell types presented decreased activity.

Specific receptor-ligand interactions also changed between control and diabetic mice. The full data set of ligand-receptor partners between clusters is available upon request. We hypothesized that changes in the proportions of cell types could be explained by shifts in genes that are involved in cell-cell interaction and expressed by other cell types. To test this, we examined all cell subset pairs for each receptor-ligand partner: whether the ligand’s expression level in one cell subset was associated with the proportions of the cell subset expressing its receptor (including autocrine communication). For instance, VEGFA upregulation by RGCs, rods, and cones during diabetes is correlated with decreased proportions of Müller glia, which express its receptor VEGFR2 that is encoded by kdr (Fig. 3B) (Spearman’s ρ = –0.26, –0.77, and –0.20, respectively). Additionally, the frequency of Müller glia, which express ephb1, was correlated with the expression by RGCs of efnb3 (Fig. 3B) (Spearman’s ρ = 0.99).

In contrast to bulk tissue gene expression, we used scRNA-seq for gene expression analysis, with the expectation that genomic information in individual cell types would help elucidate the mechanisms of DR pathology that might be hidden in bulk analysis. To identify cell-specific genes that may contribute to DR pathogenesis, we determined the top 25 up- and downregulated genes between diabetic and control retinal cells in each cell type (Fig. 4, Supplementary Fig. 5). There were several DEGs identified as multi-cell type (Fig. 5A). In comparison with their levels in control retinas, fos, madd, and pttg1 were upregulated in multiple cell types. Notably, the expression of RLBP1, which encodes CRALBP, was reduced in Müller glia but increased in several other retinal cell types in diabetic mice in comparison with that in control mice. In addition to these relatively universal DEGs, cell-type-specific DEGs could serve as selective therapeutic targets to alleviate or normalize specific pathological abnormalities of DR (Fig. 5B). For example, calcium voltage-gated channel subunit α1 H, encoded by CACNA1H, was downregulated in cone bipolar cells during DR. DEGs were used for analysis of functional enrichment with KEGG. The top 30 most significantly enriched KEGG pathways were exhibited in different cell types (Supplementary Fig. 6).

Targeting pan-retinal DEGs may offer more comprehensive and stronger therapeutic effects by normalizing multiple retinal cell types. Several DEGs were identified in multi-cell type affected by diabetes, including RLBP1, which has been widely studied in the context of autosomal recessive retinal diseases. Mutations in human RLBP1 can cause retinitis punctata albescens, Newfoundland rod-cone dystrophy, autosomal recessive retinitis pigmentosa, and Bothnia dystrophy (18). However, little is known regarding its role in DR, possibly because previous bulk measurements could not distinguish changes in RLBP1, which is downregulated in DR Müller glia but upregulated in other retinal cell types.

We hypothesized that upregulated RLBP1 expression in the Müller glia could serve as a target to mitigate neurovascular degeneration in DR. To our knowledge, we first confirmed the significant downregulation of CRALBP in diabetic retinal tissue using immunofluorescence staining and Western blotting (Fig. 6A and B). After verifying that intravitreal AAV injection efficiently delivered the transgenes to the Müller glia (via the Müller glia–specific ShH10Y.GFP [15] [Fig. 6C]), we used the AAV vectors to express RLBP1 specifically in the Müller glia of db/db mice and recorded the fasting body weights and blood glucose levels of all four groups (Supplementary Table 2). We also examined whether RLBP1 overexpression could affect glial dysfunction and proinflammatory cytokine induction. Glial responses were assessed with the gliotic marker GFAP. As shown in Fig. 6D, GFAP expression was limited in astrocytes and in the end feet of Müller glia in the retinas of db/m (control) mice. However, GFAP expression in untreated db/db mice was detected in the cell processes of the Müller glia, indicating their activation or damage (Fig. 6D). However, db/db mice intravitreally injected with ShH10Y.RLBP1 displayed decreased upregulation of GFAP in the astrocytes and Müller glia (Fig. 6D). Moreover, diabetes led to increased expression of IL-1β (P < 0.001 [Fig. 6E]) and TNF-α (P < 0.001 [Fig. 6E]), whereas overexpressing RLBP1 resulted in downregulating these cytokines among db/db mice (P < 0.01 and 0.05 [Fig. 6E]). Furthermore, the number of TUNEL-positive cells, indicating DNA strand breaks, was significantly increased in the retinas of untreated db/db mice compared with that in the retinas of db/m mice (Fig. 6F). In db/db mice injected with ShH10Y.RLBP1, the number of TUNEL-positive cells in retinas significantly decreased compared with that in db/db mice administered PBS (Fig. 6F). In evaluation of neuronal dysfunction, RGC-specific marker Brn3 revealed significant attenuation in RGCs of untreated db/db mice in comparison with db/m mice (Fig. 6H). However, db/db mice intravitreally injected with ShH10Y.RLBP1 displayed increased RGCs in comparison with db/db mice intravitreally injected with PBS (Fig. 6H). In addition, the amplitudes of both a- and b-waves increased in a light intensity–dependent manner in the db/m mice. However, a- and b-waves were significantly decreased in the db/db mice that were intravitreally injected with PBS in comparison with the db/m mice from the second data point to the fourth data point. Compared with treatment with PBS, intravitreal injection of ShH10Y.RLBP1 significantly reduced attenuation of a- and b-wave amplitudes (Fig. 6I). At the 2nd, 3rd, and 4th a-wave data points, there were significant increases in the ShH10Y.RLBP1-treated group compared with the PBS-treated group. Moreover, at the 3rd and 4th b-wave data point, there was a significant increase in the ShH10Y.RLBP1-treated group in comparison with the PBS-treated group. We next assessed the impact of RLBP1 on DR vasculopathy by examining acellular capillaries (Fig. 6J). As expected, db/db mice that were administered PBS had significantly increased numbers of retinal acellular capillaries in comparison with db/m mice (P < 0.001 [Fig. 6K]), whereas db/db mice that were intravitreally injected ShH10Y.RLBP1 had a significantly decreased number of retinal acellular capillaries in comparison with those injected with PBS (P < 0.05 [Fig. 6K]). Evans blue was used to determine the extent of blood-retinal barrier (BRB) permeability. BRB breakdown and increased vascular permeability were observed in db/db mice that were intravitreally injected with PBS relative to db/m mice (P < 0.001 [Fig. 6L]). db/db mice with RLBP1 overexpression exhibited a significant decrease in vascular permeability in comparison with db/db mice that were intravitreally injected with PBS (P < 0.001 [Fig. 6L]). However, there were no significant differences in the indicators between db/db mice injected with PBS and ShH10Y.GFP (Fig. 6D–L).

Here, we described high-throughput parallel scRNA-seq analysis of retinal cells from control and diabetic mice. This enabled the characterization of genes as cell-specific markers within the complex retinal cytoarchitecture and identified genes whose expression is significantly changed in diabetes, thus providing a rich resource for the study of DR pathogenesis. Viewed holistically, our analyses of changes in cell proportions, cell-cell communication, and gene expression indicate that most retinal cell types demonstrate various degrees of sensitivity to DR.

GWAS have identified genomic loci associated with various human phenotypes; however, the specific cellular mechanisms behind these disease associations remain largely unknown. Exploring the cell-type-specific expression of DR-associated risk loci can provide deeper understanding of DR pathogenesis. Müller glia cluster 6, as well as all clusters containing rods, rod bipolar cells, cones, and vascular endothelial cells, was the most intimate relative in PDR, indicating that they may play greater roles in the disease process.

With in-depth pathological studies and improved retinal imaging to conceptualize DR, it has been found that the retinal dysfunction induced by diabetes may be viewed as an impairment in the retinal neurovascular unit. The retinal neurovascular unit contains the physical and biochemical relationship among neural cell types, glial cells, professional immune cells, and vascular cells (19). This prompted us to consider DR as a whole, with a focus on cell-cell interactions. In comparison with db/m mice, microglial activity in retinal cell-cell communication was increased in db/db mice, probably due to their diabetes-induced activation and consequent migration from the inner retinal layers into the outer plexiform and photoreceptor layers (20). However, other cell types showed decreased activity, possibly because of impaired morphology and function. Importantly, retinal neurodegeneration precedes microvascular changes in DR (21). As the Müller glia are the only cells that span the entire retina, with intimate contacts with other retinal cell types, and VEGFR2-mediated Müller cell survival is vital for retinal neuron viability in diabetes (22), we studied the association between the proportion of Müller cells expressing VEGFR2 and neurons expressing VEGF. As shown in our results, the proportion of Müller glia that expresses VEGFR2 decreased with the VEGF upregulation of RGCs, rods, and cones during diabetes. Thus, we speculate that increased VEGF expression in neurons leads to decreased Müller glia expressing VEGFR2. We also observed some ligand-receptor interactions that were altered in diabetes that may improve our understanding of DR pathogenesis. For instance, signaling between erythropoietin-producing human hepatocellular (Eph) receptors and Eph family receptor–interacting (ephrins) proteins was active in retinal tissue both among several cell types and in autocrine loops in single cell types. The Eph receptor/ephrin system controls cell fate and is typically mediated via contact-dependent intercellular communication. It has a number of roles in both development and adulthood (23,24). In retinal development, Eph receptors and ephrins participate in retinotectal mapping. They also contribute to retinal neovascularization via Eph receptor A2/ephrin A1 signaling in endothelial cells and via RGC apoptosis via Eph receptor B/ephrin B signaling in Müller glia (25,26). Based on our result, the lesser the expression of ephrin B3 (EFNB3) in RGCs, the lower the proportions of Müller glia that express EPH receptor B1 (EPHB1). This may indicate a cascade between RGCs and Müller glia via Eph receptor B/ephrin B signaling. Other Eph receptors/ephrins have not been extensively researched in retinopathy, and our study reveals their new roles in DR. Eph receptor A4 (EPHA4) affects ischemic brain injury by controlling astrocyte glutamate uptake capacity, preventing glutamate excitotoxicity, mediating neuroinflammation and tissue damage in traumatic brain injury, and promoting glial scarring after ischemic stroke (27–29). Analogously, glutamate excitotoxicity, reactive gliosis, and neuroinflammation are characteristic retinal neurodegeneration phenotypes of DR (30). Thus, the role of EPHA4 in DR should be examined further.

Single-cell data were essential for prioritization of the numerous promising therapeutic targets in the diverse cell types. There are some cell-type-specific DEGs, such as purinergic receptor P2Y12 (P2RY12). P2RY12, which was decreased in retina during diabetes in this study, is a marker in homeostatic microglia, which was lost in active and slowly expanding lesions (31). In central nervous system, microglia are resident immune cells, including the retina, and are vital for the homeostatic maintenance of the nervous and retinal microenvironment (32). Recently, microglia were proven to be essential for initiating the infiltration of immune cells in experimental autoimmune uveoretinitis (EAU) murine model (33). Microglia-Müller glia or retinal pigment epithelial (RPE) cells cross talk drives inflammation and destroy BRB in DR (34,35). Hence, it seems that future studies on P2RY12 in microglia may be valuable. Additionally, connective tissue growth factor (CTGF) is a major contributor to fibrotic responses in DR. CTGF was only significantly increased in Müller glia in our study, suggesting that these cells are the major source of CTGF in DR. DEG enrichment analysis also revealed significant pathways and gene sets that changed in different cell types with diabetes. For instance, mitophagy-related genes were strongly enriched in RGCs, consistent with results from studies on cell differentiation and glaucoma (36–38). This indicates that further research into RGC mitophagy is needed to shed light on DR pathogenesis and provide novel treatment strategies.

Pan-retinal DEGs were identified in multi-cell types. Notably, RLBP1 one of the most robust DEGs across retinal cell types, exhibited opposing trends between the Müller glia and other retinal cell types. Müller glia form the supporting architectural structure across the entire thickness of the retina and can be viewed as the core of the retinal neurovascular unit, as they provide an anatomical connection between retinal blood vessels and neurons (19). The three primary functions of Müller glia are 1) uptake and recycling of retinoic acid compounds, neurotransmitters, and ions; 2) blood flow regulation and BRB maintenance; and 3) control of retinal metabolism and nutrients. Müller glia play a noteworthy role in DR, causing chronic inflammation, apoptosis, vascular leakage, and neovascularization (39). After selectively overexpressing RLBP1 in diabetic Müller glia, we observed decreased gliosis and protected retinal capillaries and neurons. RLBP1 is an 11-cis-retinoid-binding protein that is abundantly expressed in both RPE cells and Müller glia and is vital for the retinal visual cycle. Phototransduction is triggered when the 11-cis-retinal chromophore of a photopigment molecule in the outer segment of a vertebrate rod or cone captures a photon and is isomerized to its all-trans form. Regeneration of light detection requires reconverting the “bleached” all-trans retinal back into a light-sensitive 11-cis retinal (40). In RPE cells, RLBP1 contributes 11-cis-retinal to rods and cones, while cones also depend on an 11-cis chromophore present in Müller glia, which explains the characteristic rapid dark adaptation and saturation resistance of cones (41). Accumulating evidence suggests that visual dysfunction and retinal neurodegeneration precede retinal microvascular degeneration (21). Impaired retinoid metabolism during diabetes is characterized by 11-cis-retinal deficiency, resulting in increased unbound opsin that constitutively activates the phototransduction cascade, causing photoreceptor damage and exaggerating oxidative stress (42). As 11-cis-retinal is highly unstable, Malechka et al. (42) showed that the administration of the relatively stable 11-cis-retinal isomer 9-cis-retinal to a model of type 1 diabetes rescued visual pigment formation and decreased oxidative stress and retinal neuron apoptosis. How RLBP1 overexpression in Müller glia alleviates diabetes-driven retinal neurovascular degeneration remains a matter of speculation; however, it is possible that the normalization of retinoid metabolism results in decreased oxidation stress and apoptosis, while maintaining Müller glia homeostasis to benefit the entire retina.

Our study has limitations. Because of the differences between mouse and human retinas, mouse models only partially recapitulate the pathological complexity of patients with DR, and further research will be necessary to translate this information in humans. In addition, retinal scRNA-seq at different time points of db/db mice should be added for full understanding of the occurrence and development of the disease. Another challenge is the accurate detection of rare cell types using scRNA-seq. Distinguishing rare cell types from massive scRNA-seq data sets and building their precise signatures are challenging tasks (43), and estimating rare cell proportions is particularly difficult with small sample sizes because of increased stochasticity. We hope to further improve on the detection accuracy of rare cell types in future studies.

To summarize, this study identified DR-associated genes, cell-cell communication pathways, and retinal cell types, laying the foundation for further DR pathogenesis studies and for the development of novel treatment strategies for DR and related disorders.

Acknowledgments. The authors thank OE Biotech Co., Ltd (Shanghai, China) for providing single-cell RNA-seq.

Funding. This work was supported by grants from the National Key R&D Program of China (grant 2016YFC0904800 and 2019YFC0840607), National Science and Technology Major Project of China (grant 2017ZX09304010), National Natural Science Foundation of China (81870667 and 81800799), Shanghai Medical Excellent Discipline Leader Program (2017BR056), and Shanghai Municipal Education Commission–Gaofeng Clinical Medicine Grant Support Program (20161426).

The funding organizations had no roles in the design, conduct, analysis, or publication of this study.

Duality of Interest. No potential conflicts of interest relevant to this article were reported.

Author Contributions. T.N. and K.L. designed the research. T.N., X.S., M.Z., and Y.W. performed the research. J.F., X.X., and Y.W. analyzed the data. T.N., S.Z., and K.L. wrote and revised the manuscript. K.L. obtained funding and supervised the study, and all authors read and approved the final manuscript. K.K. is the guarantor of this work and, as such, had full access to all the data in the study and takes responsibility for the integrity of the data and the accuracy of the data analysis.