Metagenome Investigation of Ocular Microbiota of Cataract Patients With and Without Type 2 Diabetes

Abstract

Purpose: Our objective was to investigate differences in the ocular surface bacterial composition in cataract patients with and without type 2 diabetes (T2D).

Methods: Twenty-four diabetic patients with cataracts (group D) and 14 sex- and age-matched patients with age-related cataracts (group N) were recruited for this study. All samples underwent DNA extraction, fragmentation, purification, library construction, and metagenomic sequencing.

Results: The overall conjunctival sac bacterial composition was similar between group D and group N, as determined by alpha diversity and beta diversity. Nevertheless, significant differences were observed in the relative abundance of specific bacteria. At the phylum level, group D had a significantly lower abundance of Chlamydiae, Tenericutes, Chloroflexi, Cyanobacteria, Cossaviricota, Chytridiomycota, Artverviricota, Zoopagomycota, Peploviricota, Deinococcus-Thermus, Preplasmiviricota, and Nucleocytoviricota. At the genus level, group D had a significantly lower abundance of Chlamydia, Mycoplasma, Salmonella, Chryseobacterium, Roseovarius, Desulfococcus, Kangiella, Anaerococcus, and Idiomarina but a significantly higher abundance of Parabacteroides, Phocaeicola, and Sphingomonas. Bacteria such as Aquificae, Parabacteroides, Flavobacterium, Austwickia, Aquifex, Tenacibaculum, and Chryseobacterium in group D and Tenericutes, Chlamydiae, Porphyromonas, Mycoplasma, Chlamydia, Kangiella, Idiomarina, Roseovarius, Aliiroseovarius, and Desulfococcus in group N could be used as conjunctival sac biomarkers, according to the linear discriminant analysis effect size. Gene Ontology functional annotation indicated that bacterial catalytic activity, metabolic processes, locomotion, virion, and reproduction were enriched in group D, while immune system processes were enriched in group N. In addition, the top 30 differentially expressed virulence factors (VFs) were all more enriched in group D.

Conclusions: The bacterial composition was similar between the two groups. Several bacterial strains that were reported beneficial in gut were decreased, and pathogenic bacteria were increased in T2D. Furthermore, group D had more active bacterial terms and increased VF expression, suggesting that the susceptibility of diabetic patients to infection is closely related to functional changes in the ocular surface flora. Our conjunctival microbiota atlas provides a reference for investigating ocular complications related to diabetes.

Translational Relevance: The altered composition and functional profile of the ocular microbial community in diabetic patients offer evidence for the need to prevent infection during cataract surgery.

Introduction

The National Institutes of Health initiated the Human Microbiome Project in 2007 with the aim of identifying the microorganisms present in the human body and understanding how these change in relation to disease. Conjunctiva, as an ocular tissue that is directly exposed to the environment, has a complex microbiota that includes pathogenic microorganisms, opportunistic pathogens, and symbiotic bacteria.^1,2 We previously identified the core microbiome of conjunctival sac microbes in healthy individuals.³ There have been numerous reports linking the commensal microbiota in the conjunctival sac to ocular diseases such as conjunctivitis, keratitis, uveitis, dry eye, and diabetic retinopathy,^4–9 demonstrating the interaction between the host ocular disease status and the microbial community on the ocular surface. Clinical practice has shown that managing the ocular surface microbiota before, during, and after internal eye surgery can reduce or prevent postoperative complications, such as endophthalmitis.^10,11 Hence, it is important to examine the human ocular microbiota for better clinical prevention and treatment of eye diseases. 

According to the International Diabetes Federation, the number of people worldwide with diabetes mellitus was estimated to be 537 million in 2021 and is projected to reach 700 million by 2045.^12,13 Furthermore, there is a growing population of younger age with type 2 diabetes (T2D).¹⁴ Diabetes has a significant impact on the risk of infection.^15–18 Hyperglycemia could compromise the host immunity to provide a favorable environment for bacterial growth, both on the ocular surface and in tears, and disrupt the community of native bacteria.^2,19,20 

Metagenomics has been widely applied with advancements in sequencing technology for examining changes in the commensal microbiota in disease conditions.^21,22 This approach enables determination of both the taxonomic composition and functional capabilities of microbial communities, as well as the interaction between the host and the microbiome. A previous study utilizing 16S ribosomal RNA (rRNA) sequencing found notable disparities between the ocular surface microbiota of diabetic patients and healthy individuals.²⁰ Our study went beyond simply determining bacterial composition, offering insights into the functional changes of the microbial community in diabetic patients. The results shed light on the direct impact of commensal bacteria on the risk of diabetes-related infections during cataract surgery. 

Materials and Methods

Ethics Statement

This study was reviewed and approved by the Ethics Committee of Eye Institute of Shandong First Medical University (2019-30). Informed consent was obtained from all the participants. Investigations and measurements were performed in accordance with the Declaration of Helsinki. 

Inclusion and Exclusion Criteria

Patients were diagnosed with T2D when they satisfied any of the 2017 US American Diabetes Association diagnostic criteria: (1) glycosylated hemoglobin ≥6.5%, (2) fasting blood glucose ≥7.0 mmol/L, and (3) random blood glucose with typical hyperglycemia or hyperglycemia crisis ≥11.1 mmol/L. Fasting was defined as no caloric intake for at least 8 hours. 

A total of 38 patients undergoing cataract surgery were recruited at the Qingdao Eye Hospital. Patients with high intraocular pressure, glaucoma, and lens dislocation were exempted from the survey. Exclusion criteria included the administration of antibiotics, corticosteroids, or nonsteroidal anti-inflammatory drugs by oral or eye drops, concurrent corneal contact lenses, ocular surgery history, and ocular or systemic diseases that might interfere with the ocular surface indigenous microbiota. In our study, diabetic patients did not exhibit hyperglycemia associated with type 1 diabetes mellitus (DM) or other types of DM. 

Sample Collection

Samples were collected following a standardized protocol in an ophthalmic treatment room sterilized by ultraviolet irradiation. Conjunctival swab samples were obtained from one randomly selected eye or both eyes from patients scheduled for cataract surgery. After topical anesthesia with 0.4% oxybuprocaine hydrochloride eye drops (Santen, Osaka, Japan), conjunctival swab specimens were collected from the upper and lower conjunctival sac caruncle and fornix conjunctiva by gently rubbing twice with a sterile disposable swab (Guangzhou Weimi Bio-Tech, Guangzhou, China). Each swab was immediately placed into a sterile tube and stored at –80°C before DNA extraction. Strict sterile procedures for sample collection and transport were followed to prevent contamination. Sterile swabs subjected to the same sampling conditions were used for DNA extraction to monitor DNA contamination. 

DNA Extraction and Metagenome Sequencing

DNA was extracted from the conjunctival swab with a DNeasy PowerSoil Kit (Mobio, Carlsbad, CA, USA), and the DNA concentration and purity were measured with NanoDrop One (Thermo Fisher Scientific, Waltham, MA, USA). DNA fragmentation was performed using S220 Focused-ultrasonicators (Covaris, Woburn, MA, USA), and purification was performed by Agencourt AMPure XP beads (Beckman Coulter Co., Brea, CA, USA). Library construction was performed utilizing the TruSeq Nano DNA LT Sample Preparation Kit (Illumina, San Diego, CA, USA). Polymerase chain reaction cycles used were 16 for library construction. 

Bioinformatics Analysis

Libraries were sequenced on the Illumina Novaseq 6000 sequencing platform, and 150-bp paired-ended reads were generated. Raw reads were trimmed and filtered with Trimmomatic (v0.36).²³ Reads were compared to the host genome using Bowtie2 (v2.2.9)²⁴ to remove the host contamination. Metagenome splicing was assembled using MEGAHIT (v1.1.2) after obtaining valid reads.²⁵ Contigs with a length of 200 bp or 500 bp were formed (Scaftigs) and retained. Prodigal (v2.6.3)²⁶ was used to predict open reading frames using assembled scaffolds and translation into amino acid sequences. Non-redundant (NR) gene sets were constructed for all predicted genes using CDHIT (v4.5.7).²⁷ The clustering parameters were 95% identity and 90% coverage for the sequence. In each cluster, the longest gene was selected as the representative sequence of the gene set. After obtaining the gene set's representative sequence, each sample's clean reads were compared with the NR gene set separately (95% identity) using bowtie2 (v 2.2.9) to count the abundance information of the gene in the corresponding sample. Taxonomy was derived from the taxonomy database of the National Research Library. Abundance statistics were applied at each taxonomic level (domain, kingdom, phylum, class, order, family, genus, and species) to construct the abundance profile. Annotations were performed with the NR, Kyoto Encyclopedia of Genes and Genomes (KEGG), The Clusters of Orthologous Groups (COGs), SWISS-PROT, and Gene Ontology (GO) databases, with an e-value of 1e-5 using DIAMOND (v 0.9.7). Hmmscan (v3.1b2) was used to compare gene sets with the CAZy database to obtain information on carbohydrate-active enzymes corresponding to genes. Then, carbohydrate activity was calculated as the sum of gene abundances according to carbohydrate-active enzyme abundance. R package (v 3.2.0; Indianapolis, IN, USA) was used for principal component analysis (PCA) analysis, species or functional abundance spectra plotting, principal coordinates analysis (PCoA) distance matrix calculation, and graphical analysis. The difference between the two groups was compared by the Wilcoxon test. Linear discriminant analysis effect size (LEfSe) was applied to analyze variance in species abundance or functional abundance spectra. Metagenomic sequencing and bioinformatics analysis were performed by OE Biotech Co., Ltd. (Shanghai, China). 

Statistical Analysis

R software (v 3.2.0), GraphPad Prism 5.0 (GraphPad Software, San Diego, CA, USA), and SPSS 20.0 software (SPSS, Inc., Chicago, IL, USA) were used for statistical analyses. Wilcoxon and Mann–Whitney U test were used to compare the difference between the two groups. P < 0.05 was considered statistically significant. 

Results

Participant Characteristics

This study included 38 participants with an average age of 67.89 ± 7.48 years (mean ± SD). Out of these, 24 were diagnosed with cataracts and T2D (group D), consisting of 12 males and 12 females. The control group (group N) consisted of 14 age-related cataract patients, including 7 males and 7 females. Table 1 displays the demographic information of the two groups. No significant difference in age was observed between the two groups (P > 0.05, Mann–Whitney U test). 

Table 1.

 

Demographic and Clinical Characteristics

 

Overview of the Microbiome

Two sterile swabs subjected to the same sampling conditions were prepared as negative control. No visible DNA band was generated from these negative control samples, indicating that the sampling process was effectively controlled and there was little possibility of environmental contamination in metagenomic sequencing (Supplementary Fig. S1). Host contamination was eliminated prior to library construction, and the percentage of host reads is presented in Supplementary Table S1. For samples in group D and group N, an average of 31.62 phyla and 306.87 genera were identified from each subject. Only 0.47% of the taxa identified in group N were fungi and 1.60% were viruses. For group D, 0.39% were fungi and 1.45% were viruses. Over 90% of the taxa identified in both groups were bacteria (Supplementary Fig. S2). Therefore, we focused on analyzing bacterial spectrum in this study. 

Bacterial Compositions Were Similar in Bilateral Conjunctival Sacs in Group N

In our study, we followed a two-phase approach (Fig. 1A). During the first phase, we enrolled four patients with age-related cataracts and sequenced both of their eyes (N1–N8) to check for any disparities in bacterial composition between the two eyes of the same person. N1 and N2, N3 and N4, N5 and N6, and N7 and N8 correspond to the two eyes of the same individual, with N1, N3, N5, and N7 representing the left eye group (group L) and N2, N4, N6, and N8 representing the right eye group (group R). Similar to the study by Cavuoto et al.,²⁸ we found no significant difference between the two eyes of the same individual at both phylum and genus levels (P > 0.05, t-test) (Figs. 1B, 1C). Therefore, we assumed that there were similar bilateral conjunctival sac commensal bacteria in our cataract cohort. In the second phase, we only sequenced one randomly selected eye from each individual. 

Figure 1.

 

(A) The experimental design of our study. (B) The composition and relative abundances of the top 15 phyla in N1 to N8. (C) The composition and relative abundances of the top 15 genera in N1 to N8. The relative abundance was presented in the stack.

Difference of Ocular Surface Commensal Bacterial Communities Between T2D and Nondiabetic Cataract Patients

To identify the metagenomic signatures of the ocular surface related to T2D, we compared the bacterial spectrum and dominance between group N and group D. Actinobacteria was the most abundant phylum in both groups (Fig. 2A). The top six phyla shared by these two groups were the same, being Actinobacteria (average abundance of 59.48% for group N and 66.17% for group D), Proteobacteria (14.98% and 13.87%), Firmicutes (10.31% and 8.50%), Bacteroidetes (4.59% and 2.80%), Kitrinoviricota (0.73% and 1.12%), and Chlamydiae (0.85% and 0.30%) (Fig. 2A). The top 12 differentially expressed phyla between the two groups were Chlamydiae, Tenericutes, Chloroflexi, Cyanobacteria, Cossaviricota, Chytridiomycota, Artverviricota, Zoopagomycota, Peploviricota, Deinococcus-Thermus, Preplasmiviricota, and Nucleocytoviricota. These phyla were found in significantly lower amounts in group D (P < 0.05, Wilcoxon signed-rank test) (Fig. 2B). 

Figure 2.

 

Differences in the microbiota composition between group D and group N. (A) The composition and relative abundances of the main phyla in the two groups. (B) The difference in the relative abundance of the top 12 phyla between group D (orange) and group N (blue). (C) The composition and relative abundances of the main genera in the two groups. (D) The difference in the relative abundance of the top 12 genera between the two groups. *P < 0.05. **P < 0.01. ***P < 0.001.

At the genus level, Corynebacterium was the most abundant genus in both groups (Fig. 2C). The top 10 genera shared by both groups were the same, being Corynebacterium (34.32% for group N and 44.50% for group D), Cutibacterium (12.75% and 8.86%), Escherichia (6.03% and 3.73%), Mycobacteroides (4.95% and 3.47%), Staphylococcus (3.43% and 3.22%), Francisella (2.09% and 2.61%), Propionibacterium (2.88% and 1.99%), Lawsonella (1.71% and 2.56), Bacillus (2.71% and 1.82%), and Streptomyces (0.26% and 2.41%) (Fig. 2C). Group D had a lower abundance of Chlamydia, Mycoplasma, Salmonella, Chryseobacterium, Roseovarius, Desulfococcus, Kangiella, Anaerococcus, and Idiomarina but a high abundance of Phocaeicola, Parabacteroides, and Sphingomonas (P < 0.05, Wilcoxon signed-rank test) (Fig. 2D). The top 15 species within each sample are shown as histograms in Supplementary Figure S3. 

There was no significant difference in the alpha diversity of the ocular surface bacterial flora between groups N and D, as indicated by the Chao1 index (P = 0.22) (Fig. 3A) and the Shannon index (P = 0.38) (Fig. 3B). The good coverage index showed that the sequencing results were representative of the actual microorganisms and there was no statistical difference between the two groups (P = 0.29) (Fig. 3C). 

Figure 3.

 

The alpha diversity and beta diversity between group D and group N. Scatterplots of alpha diversity indices of samples from the two groups, with (A) Chao 1 index and (B) Shannon indices indicating the species richness. There was no difference between the two groups. (C) The coverage index showed that the sequencing results were representative of the actual microorganisms in the samples, and there was no statistical difference between the two groups. (D) Two-dimensional PCoA plots were constructed for beta diversity. Weighted Unifrac algorithm PCoA between group D (blue dots) and group N (red dots) at the phylum level (P = 0.715). (E) Weighted Unifrac algorithm PCoA between group D (blue dots) and group N (red dots) at the genus level (P = 0.448). There was no difference in the diversity and richness of the ocular microbiota of the two groups.

The beta diversity between the two groups was assessed by the PCoA. The samples were divided into two categories, those with similar community structures and those with distinct community structures (Figs. 3D, 3E). The PCoA plots showed that there was no significant difference in the composition of bacterial communities between groups D and N at both the phylum level (Weighted Unifrac algorithm, P = 0.715) (Fig. 3D) and the genus level (P = 0.448) (Fig. 3E). 

A comparative analysis of group D and group N microbial communities at the phylum and genus levels was performed using the LEfSe method. The log 10 linear discriminant analysis (LDA) scores for the most prevalent taxa in group N were presented on a positive scale, while the scores for the most dominant taxa in group D were presented on a negative scale. A total of 27 differentially abundant taxonomic clades were found (LDA score > 2, P < 0.01) (Fig. 4 and Supplementary Fig. S4). Aquificae was the biomarker phyla in group D and Tenericutes and Chlamydiae in group N. The unique genera in group D were Parabacteroides, Flavobacterium, Austwickia, Aquifex, Tenacibaculum, and Chryseobacterium, while in group N, they were Porphyromonas, Mycoplasma, Chlamydia, Roseovarius, Kangiella, Idiomarina, Aliiroseovarius, and Desulfococcus. 

Figure 4.

 

LEfSe analysis of group D and group N. Cladogram of the conjunctival bacterial taxa in group D (red) and group N (green). Nodes indicated taxa at different levels. The diameter of connections represented the relative abundances of taxa. p, c, o, f, g, and s are equivalent to phylum, class, order, family, genus, and species, respectively.

Gene Function Annotations of the Conjunctival Sac Microbiota

A total of 57,778 unigenes were obtained from all samples and categorized into 55 GO functional classifications (Fig. 5), including 22 biological process terms, 19 cell components terms, and 14 molecular function terms. The GO terms that were significantly upregulated in group D included catalytic activity, metabolic process, antioxidant activity, locomotion, and virion, while the immune system process was highly upregulated in group N (Table 2). These findings suggest that bacterial activities on the ocular surface were more active in group D, and the immune activity was activated in group N, which may explain why patients with T2D are more susceptible to infections. 

Figure 5.

 

GO function classification in group D and group N. A total of 55 GO functional classifications were obtained, including 22 biological process terms, 19 cell component terms, and 14 molecular function terms.

Table 2.

 

Numbers of Classification in the Top GO Terms Shared by Group D and Group N

 

For the KEGG analysis, a total of 48,725 unigenes were assigned to 290 KEGG pathways. Pathways related to aging, carbohydrate metabolism, environmental adaptation, and membrane transport were significantly enriched in group D, while pathways related to development, the digestive system, the immune system, and signaling molecules and interaction were significantly enriched in group N (P < 0.05, Wilcoxon signed-rank test) (Fig. 6A). Differentially abundant pathway clades (LDA score >2, P < 0.05) were observed (Fig. 6B), which is consistent with the GO results (Fig. 5). The highly enriched KEGG classifications for all samples included carbohydrate metabolism, amino acid biosynthesis, and energy metabolism (Fig. 6C). 

Figure 6.

 

Overview of KEGG analysis. (A) Comparison of the abundance of the top KEGG pathways between group D (red) and group N (blue). (B) The score of the linear discriminant analysis of biomarkers in group D (red) and group N (green) (LDA score >2, P < 0.05). (C) KEGG functional predictions and classifications of all identified unigenes in the two groups.

Virulence Factors Were More Enriched in T2D

The top 30 virulence factors (VFs) in both groups being displayed in Figure 7. All VFs were more abundant in group D (Figs. 7A, 7B). The top three VFs were identical in the two groups, being capsule, trehalose-recycling ABC transporter and biosynthesis and transport of phthiocerol dimycocerosate (PDIM) and phenolic glycolipid (PGL) (Figs. 7A, 7B). 

Figure 7.

 

Overview of VF in group D and group N. (A) Top 30 VF classification in group D. (B) Top 30 VF classification in group N.

Discussion

The ocular surface is exposed to the environment and hosts a group of commensal bacteria. There is growing evidence that ocular surface commensal bacteria can interact with the host immune system and consequently influence the development and prognosis of diseases.^2,29,30 Previously, we found that the bilateral microbiota of patients with fungal keratitis differed significantly, indicating alterations in flora structure may contribute to the susceptibility to fungal infection.⁷ This highlights the importance of thoroughly examining ocular surface microorganisms to gain a better understanding of the pathogenesis of ocular diseases from a new viewpoint. 

T2D is a common chronic disease with various complications.³¹ Compared to healthy individuals, patients with T2D have a higher risk of developing infections.^15,19,32 The increased susceptibility may be related to the immune system dysfunction that results from hyperglycemia.¹⁹ However, the impact of hyperglycemia on commensal bacteria is still not well understood. Studies have shown that the accumulation of advanced glycation end products in the eye due to poor glucose utilization in patients with T2D can create a favorable environment for microbial growth.^20,33 Additionally, some microorganisms become more virulent in high-glucose environments and are more likely to adhere to diabetic cells compared to nondiabetic cells.¹⁷ These findings provide evidence for the alteration of ocular surface flora structure and increased susceptibility to infections in patients with T2D. 

Our study applied metagenomics to analyze the microbiome of the conjunctival sac, offering insights into both the structural and functional changes of the microbial community in diabetic individuals. One significant observation was that the microbial composition in the eyes of cataract patients was nearly identical, as both eyes are exposed to similar external conditions. This is also consistent with findings by Cavuoto.³⁴ Thus, we randomly chose one eye per person for further analysis. Our results revealed that Corynebacterium was the most predominant genus in both groups, which aligns with our prior research.¹ 

In the present study, we found that group D had a significantly lower abundance of Chlamydia, Mycoplasma, Salmonella, Chryseobacterium, Roseovarius, Desulfococcus, Kangiella, Anaerococcus, and Idiomarina, while the abundance of Parabacteroides, Phocaeicola, and Sphingomonas was significantly higher. Anaerococcus is a gram-positive anaerobic bacterium commonly found in the human microbiota, known to positively impact host immunity and metabolism.^35,36 A decrease in Anaerococcus in diabetic patients may be linked to a weakened immune system. Phocaeicola is found in the human intestine and produces carbohydrate-active enzymes,³⁷ while Sphingomonas is a gram-negative bacterium with a cell membrane composed of sphingolipids, providing protection from antimicrobial substances.³⁸ The decrease of potentially beneficial bacteria and increase in pathogenic bacteria on the ocular surface may contribute to increased susceptibility to infections in diabetes. 

Furthermore, we found that VFs were more abundant in the diabetic group. The top three VFs were capsule, trehalose-recycling ABC transporter, PDIM, and PGL biosynthesis and transport. Previous studies have shown that the capsule could help cell wall core galactofucan synthesis, which is essential for the growth of bacteria.³⁹ Trehalose-recycling ABC transporter facilitates substrate transport across membranes of gram-positive bacteria, indispensable for the uptake of bacterial nutrients and essential for the survival and virulence of some human pathogens.⁴⁰ The higher abundance of VF suggests that the diabetic ocular surface microenvironment is conducive to the invasion, growth, replication, and transmission of infection by ocular surface or opportunistic pathogens, thereby increasing susceptibility to ocular infections. 

Previously, we used 16S rRNA sequencing to define the core bacterial microbiome of the conjunctival sac in the age-related cataract population,¹ which was identical to control subjects in this study. The 16S rRNA sequencing detected 25 phyla and 526 genera, while metagenomic sequencing detected 58 phyla and 815 genera. Metagenomic sequencing produced a more in-depth taxonomy profile. The bacterial composition was similar between the 16S rRNA study and the metagenome study for phyla with a relative abundance of more than 1%. The major phyla of the ocular surface appear to be relatively stable. The major genera, such as Corynebacterium and Staphylococcus, were also similar, but their relative abundance varied among subjects. The variations between the 16S rRNA sequencing and our metagenomic sequencing results may be due to the differences in the sensitivity of sequencing methods and in the way of sampling. Current 16S sequencing focused on community diversity and failed to produce direct functional files, while metagenomic sequencing generated more comprehensive and detailed information on the microbiota in a specific environmental sample and provided a more accurate function impact of a microbial community and its interactions with the environment.^41–44 

In previously published studies on the ocular surface microbiomes of diabetic and nondiabetic individuals,^20,45–47 out of the top five phyla in the diabetic groups, four were consistent with our findings, suggesting the major phyla on the ocular surface are relatively consistent regardless of the sequencing method. However, there were variations in the genera identified by these studies, which may be due to deference in geographical environment, sample size, study subjects, sensitivity of sequencing methods, or the way of sampling. Nonetheless, as our current study focused on the impact of diabetes on the ocular surface bacterial community, the differences in sequencing methods would not impact our conclusions, as both the diabetic group and the control group were subjected to the same experimental protocol. 

While our study has its limitations, such as the need for larger sample sizes and in vitro confirmation of differentially enriched bacteria, it provides important insights into the impact of diabetes on the conjunctival sac microbiota. In conclusion, our research presented a comprehensive view of the conjunctival sac’s microbial community in individuals with diabetes, highlighting the structural and functional differences from those without diabetes. We also shed light on how these differences may increase the risk of infection in people with diabetes. 

Acknowledgments

Supported by the Taishan Scholar Program (ts20190983 to Y.H.), the National Natural Science Foundation of China (81970788 and 82171027 to Y.H.; 82101091 to B.N.Z.), and the Shandong Provincial Natural Science Foundation (ZR2020QH140 to B.N.Z.). 

The raw sequence data reported in this study have been deposited in the Genome Sequence Archive (https://ngdc.cncb.ac.cn/gsa) with the accession number CRA010876. 

Disclosure: Z. Shao, None; X. Shan, None; L. Jing, None; W. Wang, None; W. Li, None; Z. Ren, None; B.N. Zhang, None; Y. Huang, None