key: cord-0823694-lcl3k7bo authors: Xiang, Haitao; Zhao, Yingze; Li, Xinyang; Liu, Peipei; Wang, Longlong; Wang, Meiniang; Tian, Lei; Sun, Hai-Xi; Zhang, Wei; Xu, Ziqian; Ye, Beiwei; Yuan, Xiaoju; Wang, Pengyan; Zhang, Ning; Gong, Yuhuan; Bian, Chengrong; Wang, Zhaohai; Yu, Linxiang; Yan, Jin; Meng, Fanping; Bai, Changqing; Wang, Xiaoshan; Liu, Xiaopan; Gao, Kai; Wu, Liang; Liu, Longqi; Gu, Ying; Bi, Yuhai; Shi, Yi; Zhang, Shaogeng; Zhu, Chen; Xu, Xun; Wu, Guizhen; Gao, George F.; Yang, Naibo; Liu, William J.; Yang, Penghui title: Landscapes and dynamic diversifications of B-cell receptor repertoires in COVID-19 patients date: 2021-11-04 journal: Hum Immunol DOI: 10.1016/j.humimm.2021.10.007 sha: 3371b5cd7cf5193e2767654e8f4710c3f7978745 doc_id: 823694 cord_uid: lcl3k7bo Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has caused the pandemic of coronavirus disease 2019 (COVID-19). Great international efforts have been put into the development of prophylactic vaccines and neutralizing antibodies. However, the knowledge about the B cell immune response induced by the SARS-CoV-2 virus is still limited. Here, we report a comprehensive characterization of the dynamics of immunoglobin heavy chain (IGH) repertoire in COVID-19 patients. By using next-generation sequencing technology, we examined the temporal changes in the landscape of the patient’s immunological status and found dramatic changes in the IGH within the patient’s immune system after the onset of COVID-19 symptoms. Although different patients have distinct immune responses to SARS-CoV-2 infection, by employing clonotype overlap, lineage expansion, and clonotype network analyses, we observed a higher clonotype overlap and substantial lineage expansion of B cell clones 2-3 weeks after the onset of illness, which is of great importance to B-cell immune responses. Meanwhile, for preferences of V gene usage during SARS-CoV-2 infection, IGHV3-74 and IGHV4-34, and IGHV4-39 in COVID-19 patients were more abundant than those of healthy controls. Overall, we present an immunological resource for SARS-CoV-2 that could promote both therapeutic development as well as mechanistic research. The current outbreak of SARS-CoV-2 infection is threatening global public health [1] . The scale of the humanitarian and economic impact of the COVID-19 is driving intense efforts to develop vaccines and neutralizing antibodies (NAb) against COVID-19. Therefore, understanding the principles of the B-cell responses during SARS-CoV-2 infection is of substantial importance for anti-viral vaccine and NAb development. The B-cell receptors (BCRs) are immunoglobin molecules located on B-cell surfaces to recognize and bind foreign antigens. Upon encountering their specific antigen, B-cells become activated, proliferate, and may differentiate into produce short-lived effective antibody-secreting plasma cells or long-lived plasma cells and memory cells. At the same time, BCRs undergo a process of affinity maturation, which repeats cycles of somatic hypermutation of BCRs and subsequent clonal selection leads to increased binding affinity. Comparison of BCR sequences among individuals is of great interest because repertoires may have similar features if individuals are exposed to the same pathogen, giving rise to convergent antibodies. Antibody specificity is largely determined by the IGH gene sequence used by each B-cell [2] . The recent developments in high-throughput sequencing have made it feasible to character IGH repertoire in large numbers of samples and it is increasingly being applied to gain insights into the humoral responses in healthy individuals and a wide range of diseases. This technic has also led to advances in our understanding of how the antibody repertoire changes in response to perturbation arising from initial viral infection, viral evolution, and vaccination. BCR repertoire analysis has been applied, for example, to influenza virus [3] [4] , human immunodeficiency virus [5] , varicella-zoster virus [6] , and dengue virus [7] . However, the dynamics of antibody response elicited by SARS-CoV-2 infection remain to be determined. To understand how B-cell immune repertoire changes over time during SARS-CoV-2 infection, we obtained IGH repertoires from the peripheral blood samples which were collected multiple times from five COVID-19 patients. We classified sequences into clones by lineage clustering analysis and tracked changes in the following characteristics: unique number of CDR3, Shannon index, the number of high-frequency clones, cumulative frequency of the top 100 clones, and V, J-gene segment usage at a different course of time. We also conducted clonotype overlap, lineage expansion, and CDR3 sequence network structure to explore the similarity and varying trends of IGH repertoire status at different time points. Overall, during the high degree of heterogeneity in B-cell clonal dynamics among patients, key common patterns were observed, i.e. the higher clonotype overlap and substantial lineage expansion of COVID-19 patients after 2-3 weeks of onset of illness. platform with 400 bp single-end reads. The single-end 400 bp reads of the samples from the five COVID-19 patients sequenced at the MGISEQ-2000 platform were analyzed by IMonitor [9] software. Briefly, the raw reads were filtered and trimmed by SOAPnuke [10] to remove low-quality reads and adapter sequences. Reference germline gene segments were collected from the IMGT database [11] . BLAST [12] program was used to align all clean reads to a reference, and then a second alignment procedure was executed to improve the alignment accuracy. Subsequently, structure analysis was performed to deduce the structure of BCR molecular, including variable (V), diversity (D), and joining (J) segments usage, the nucleotide and amino acid sequence of the CDR3 region, and the deletions and insertions at rearrangement sites. Finally, several characters reflecting the status of the immune repertoire were counted, such as V-J pairing, V/J usage, CDR3 sequence frequency, and CDR3 length distribution. We define sequences belonging to the same clonal lineage as those sharing the same V and J germline gene combination, having similar CDR3 length, and at least 90% nucleotide sequence identity in the CDR3 region. According to this definition, we first grouped sequence with identical V and J germline gene origination and same CDR3 length, following clustering those groups by CD-HIT (version 4.6) [13] with the following parameters: -c 0.9, -G 1, -b 20, -d 0 and -n 9. To eliminate the impact of the difference in data volume on clonal diversity, we normalized each of the IGH clones to a unanimous number. Firstly, we counted an effective data size (with an "in-frame" tag in the IMonitor result) for each sample and removed some samples with significantly fewer effective data than others. Secondly, normalized data with the same size were extracted randomly from each sample by an in-house Perl program. There are several commonly used parameters reflecting the features of the BCR repertoire, including unique CDR3 number, Shannon index, high-frequency clone number, and cumulate frequency of the top 100 clones. we are using an in-house Perl script to count those parameters base on normalized data from each sample. The Shannon index is frequently used as a measure of repertoire diversity, therefore we calculated the Shannon index of observed clones from normalized data using the following function: where N is the total number of unique CDR3, and p(i) is the frequency of a single CDR3 in the BCR repertoire. After data normalization, gene usage frequency of V and J segments were counted and a k-means hierarchical clustering process was executed by R package "hclust" to compare V/J segment usage similarity both within and among individuals. Repertoire overlap analysis is a similarity measure to evaluate the correlation of overlapping clonotypes between samples. The Jaccard index, also known as intersection over union and the Jaccard similarity coefficient, is commonly used to gauging the similarity and diversity of different repertoires. To explore the impact of onset time on the status of the BCR repertoire, we counted the Jaccard index to assess the association between the repertoire at different onset times. The index is calculated as follows: where |A | denotes the total number of unique CDR3 in sample A, |A∩B| represents the intersection of the unique CDR3 number in samples A and B. We calculate pairwise overlap rate as Jaccard index multiplied by 100, and divided the pairwise overlap rate into 3 clusters: the first cluster is "2-3 week adjacents" including data of time points 12-13 days versus 15-16 days and 15-16 days versus 18-19 days after onset; the second cluster is "other adjacents" which contains two adjacent time points except "2-3 week adjacents"; the third cluster is "non-adjacent" which contains all of the paired non-adjacent time points. A two-sample Wilcoxon test (Wilcoxon test in R) was used to compare the differences between clusters. A lineage expansion analysis was performed based on the normalized data. Firstly, a "join samples" process was executed to join lineage information at all time points together, inside each patient, then the normalized lineage frequency was computed as the geometric mean of lineage frequencies that comprise a given joint lineage in intersected time points. if the lineage is missing. Its frequency is set to 1e-9. Only the most abundant CDR3 sequence was selected as a representative clonotype in each lineage. Secondly, an R program originally generated by VDJtools [14] with some modifications was used to display the expansion and contraction of lineages. The network is generated based on the BCR CDR3 amino acid similarity, each vertex in the network represents a CDR3 amino acid sequence, the size of the point represents the abundance of the clone, any two vertices are connected with an edge if they meet the following three criteria, including the same length, the same V and J gene usage, and only one amino substitution. The diagram was plotted by the "igraph" R package (https://igraph.org/redirect.html). The parameters of the network, including vertices, cluster number, and diameter of the largest cluster, are calculated using "V", "clusters", "diameter" functions in the "igraph" package. The Gini coefficient is used to measure the inequality in distribution, and we calculated the Gini coefficient of clone frequency distribution and cluster clone number distribution using the following function: where y(i) is the frequency of a single clone; and n is the total unique clone number in the calculation of clone frequency Gini coefficient; y(i) is the clone number of a single cluster, and n is the total cluster number in the calculation of cluster Gini coefficient. Among the five patients diagnosed with COVID-19, the 85-year-old patient developed severe symptoms, and the other four patients aged between 15 to 45 developed mild symptoms (Table 1) . Blood samples were collected on days 2-5 after diagnosis with 2-3 days intervals until they were tested negative for SARS-CoV-2 on days 17-19. All of them were discharged from the hospital upon recovery. Notably, 3 of the 5 patients are from the same family, that is, Patient 2 is the son of Patient 1 and Patient 3. During the process of COVID-19 diagnosis, several clinical indicators were detected by routine blood examination. At least 3 of the 5 patients showed pronounced lymphopenia. The other 2 have too few tested indicators and can't be reasonably counted. Along with the progression of COVID-19, the percentage of T lymphocytes went through a process of decreasing initially and then rebouncing. The phenomenon of lymphopenia in COVID-19 patients in our results is highly consistent with previous studies. It is well established that lymphopenia is a major immunological abnormality that occurs in COVID-19 patients [15, 16, 17, 18] . However, we found that the percentage of B lymphocytes showed an opposite trend to T cells: among these patients, the percentage of B cells in Patients 3, Patient 4, and Patient 5 reached the highest value until 14 th , 6 th and 17 th day after the onset of the disease, respectively (Supplementary Table 1 ). We performed longitudinal IGH repertoire sequencing to examine the global effects of SARS-CoV-2 infection on the B-cell compartment (Fig. 1a) . To enable direct comparisons among time points of each patient, datasets of all time points of each patient were sub-sampled to a unanimous number of clonotypes. We investigated the diversity of the IGH repertoires through four parameters, which are unique CDR3 number, Shannon index, high-frequency clone number (defined as frequency > 0.1%), and cumulate frequency of top 100 clones with the highest frequency. These parameters fluctuate substantially throughout time with no common tendencies among patients (Fig. 1b,e) , which may be due to the dramatic change of the predominant clones and lineages. Next, we thought to test whether there are preferences of V and J gene segments during SARS-CoV-2 infection. By using a hierarchical clustering strategy, we compared the V and J gene usage frequency of COVID-19 patients with that of healthy donors. We found that the most frequently observed V segments were IGHV3-23, IGHV3-30, IGHV4-34, IGHV4-39, and IGHV4-59 (Fig. 1f, g) , which is similar to previous studies [19, 20] . Additionally, the germline gene usage of IGHV3-74, IGHV4-34, and IGHV4-39 in COVID-19 patients were more abundant than those of in healthy controls. Preferred J gene segments showed a high degree of consistency among all time points in both patients and healthy donors (Fig. 1g) , with the most frequent J segment is IGHJ4, which is consistent with the previous report [21] . For each individual, the frequency of distinct gene segments usage over the infection time course was fairly stable, shown as high correlations in the heatmap of Fig. 1f and g, suggesting that the overall V, J segment useage is not markedly altered by SARS-CoV-2. This phenomenon may reflect that the human immune system is highly variable between individuals but stays relatively stable over time in a given person [22] . gene family (Fig. 2b) . Notably, after extracting the clones of the IGHV3 gene family separately, we found that this part of the clones had significantly higher SHM rate than other gene families, in all patients (Fig. 2c) . Neutralizing antibodies are responsible for defending cells from pathogens. Upon triggering by infections, the body can produce neutralizing antibodies as part of its immune response. To examine the potential of our IGH repertoire to provide basic support for the development of neutralizing antibodies, we collected 15 CDR3 amino acid sequences of antibodies with neutralizing potential from two recently published studies, and conducted comparison analysis with our datasets. Comparing with our datasets, , we found that Dr. Wu's [23] study had fewer but higher similar CDR3s than Dr. Cao's [24] study (Supplementary Fig 1) . Specifically, examing the sequences of antibody BD-494 and BD-501 from Dr. Cao's study, we found 14 and 1 clone, respecitively, having identical CDR3 amino acid sequence (100% identity score reported by blastp) in our dataset, (Supplementary Table 2 ). However, there is no clone in our data has exactly the same CDR3 amino acid sequence as the antibody reported in the Wu's study, the highest sequence similarity is 93.75%. Notably, we detected a total of 226 and 497 unique CDR3s amino acid sequences similar to the anti-SARS-CoV-2 single-domain antibody n3010 and n3088 isolated by Wu et al, respectively. Among them, Patient 2 showed the highest enrichment of CDR3s similar to those antibodies (Fig. 3a) . Further analysis to explore the usage patterns of germline V genes shows that clones with CDR3 amino acid sequences similar to potential neutralizing antibodies has a clear preferential usage of Ig heavy chain variable region 3 (IGHV3) and 4 (IGHV4) subfamily genes. It is worth noting that there is a branch also had a bias toward IGHV1-24 in clones similar to n3010, and a branch that is biased towards IGHV2-70 in clones similar to n3088. (Fig. 3b) . Next, we performed phylogenetic sequence analysis on two subsets of clones similar to antibody n3010 and n3088 respectively. Phylogenetic analysis revealed that the CDR3s belonging to the IGHV1 and IGHV2 gene families clustered into the unique clusters respectively, while the CDR3s of the IGHV3 and IGHV4 gene families were fused together (Fig. 3c) . To better understand the potential mechanism behind changes in IGH repertoire diversity, we examined how the composition of the clonotype changed during the infection. We evaluated the level of clonotype overlap between any two-time points of each patient by using an all-versus-all pairwise overlap analysis. In general, the overlap rate between different time points ranges from 0.36% to 6.57% (average 2.13%) in each individual ( Fig. 4a-e) . As expected, samples from adjacent sampling batches share more clones than those from other batches. A notable rule that shows in almost every patient is that the overlap rate of adjacent batches peaks at around 2-3 weeks. We found that "2-3 week adjacents" cluster shared significantly more clones than "other adjacents" (p=7.78 x 10 -3 ) and "non-adjacent" (p=4.52 x 10 -5 ). (Fig. 4f) . This indicates that, compared with other time points, the composition of the IGH immune repertoire within individuals is more similar during 2-3 weeks after the onset. The 2-3 weeks after onset is likely a critical stage of B-cell immune response. Next, we investigated the relative degree of clonal expansion. To measure the alterations in lineage frequency over time, we first defined clonotypes as a lineage if they meet 3 criteria: sharing the same VH and JH germline gene segments, having the same nucleotide length in the CDR3 region, and at least 90% nucleotide identity in the CDR3 region. Using this strategy, we computed the geometric mean of linage frequencies of all samples and performed a lineage tracking analysis using VDJtools. Here, we defined the highly expanded lineages as the 100 most abundant lineages (top 100 lineages). We observed an evident lineage expansion process over time. As shown in Fig. 5a- BCRs with similar CDR3 amino acid sequences are supposed to be gathered in clusters based on antigen-driven SHM. To investigate the degree of clone convergence during the disease, we build the BCR repertoire CDR3 sequence network for each individual, using data of different times after the onset of the disease based on a specific approach [25] . Networks of all the 5 patients after symptom onset have been plotted in Supplementary Fig. 2 and the network of the 41-year-old female Patient 1 at different time points is demonstrated in Fig. 6a since its relatively high number of CDR3 clones (vertices). Here, a clear clonal clustering phenomenon could be seen from the network of this patient on days 18-19 (Fig. 6a) . Subsequently, three network-related parameters, namely diameter of the largest cluster, Gini coefficient of cluster clone number, and Gini coefficient of clone frequency, were counted and summarized in Fig. 6b-d. There was a similar trend that first goes up and peaks around 2-3 weeks and falls afterward in these three parameters of almost all patients. Each parameter of the network has its implication, firstly, the diameter parameter exhibit in Fig. 6b is the length of the longest geodesic (edges between any two reachable vertices) of a network. The result indicates a higher evolution level of BCR clones in the repertoire from 2-3 weeks. Secondly, the Gini coefficient of cluster clone number which represents the level of clonal clustering also peaked around 2-3 weeks. Thirdly, the Gini coefficient of clone frequency which indicates the skewed expansion of several clones was very much synchronized with the former parameter, despite they were relatively independent. Moreover, the most abundant clones usually don't appear in the largest cluster (Fig. 6d, Supplementary Fig. 2 To obtain a more detailed understanding of the dynamics and sophistication of B-cell immune response in COVID-19 patients, we applied a comprehensive analysis of IGH repertoire during the disease process. Given the advantages of longitudinal data in profiling the fast repertoire dynamics through time in individuals responding to an immune challenge [19, 26] , we collected a set of longitudinal immune repertoire data from 5 COVID-19 patients and conducted a detailed analysis of their B-cell immune repertoires. The simultaneous profiling of clonotype overlap, lineage expansion, and clonotype network enabled us to find that a period of 2-3 weeks after symptom onset is of great importance to B-cell immune response. In short, these results reveal a clear dynamic landscape of BCR repertoire in COVID-19 patients, as well as the characteristics of the immune repertoire at different time points after symptom onset. V, J gene usage analysis revealed that each individual's frequency of different V and J gene segment usage stays relatively stable over time (Fig 1f) . Similar results have also been reported by Chen Wang et . al [6] suggesting that V, J segment use is not markedly altered by SARS-CoV-2 infection. However, the overall use of V components among patients shows a differentiated preference. As Petter Brodin mentioned [22] , this may reflect that the human immune system is highly variable between individuals but relatively stable over time within a given person. We observed a preference of the use of V gene segments inconsistent with the single-cell BCR repertoire research by Wen et al [27] , but consistent with Jacob's [28] bulk BCR repertoire study. For example, Wen's group reported that IGHV3-23 was over-represented in patients compared with healthy controls, while both our and Jacob's result showed that IGHV3-23 is highly represented both in patients and healthy controls. Additionally, we and Jacob's group both found that IGHV4-34 is more frequently used in COVID-19 patients compared with healthy donors. We believe that, given the advantages of large amounts of sequencing data, bulk BCR sequencing may be more suitable for BCR gene usage frequency assessment. Nevertheless, with only five patients, the power of our study is almost certainly constrained by the small sample size. Larger samples and more sequencing data need to be collected to help further understand these differences. We found from the longitudinal immune repertoire data that all patients' BCR repertoire composition in 2-3 weeks of onset is different from that in other periods. There are significantly more common clones between adjacent time points and evident lineage expansion during 2-3 weeks. Additionally, network analysis revealed that each patient's average degree and diameter indicators peaked around 2-3 weeks after onset. This timeline of repertoire change is consistent with an antigen (virus) induced B-cell differentiation and maturation. The maturation process is on par with other viruses, such as influenza viruses [29] , and a similar phenomenon and timeline has also been observed in patients with the Middle East respiratory syndrome (MERS) [30] . We speculated that the higher clonotype overlap and substantial lineage expansion of COVID-19 patients after 2-3 weeks of illness may be associated with strong neutralizing antibody responses, accompanied by a large expansion of virus-specific B-cell clones and finally causing Bcell clone convergence. Many studies have shown that the monoclonal antibodies targeting viral surface antigens have therapeutic and preventive effects on infectious diseases such as HIV, Ebola, and MERS [31, 32, 33] . Their safety and effectiveness in patients have been confirmed in many clinical trials [34, 35] . With the efforts of scientists from various countries, recent positive progress has been made in the development of anti-SARS-CoV-2 drugs [36, 35, 37, 38, 39, 23] . We compared the IGH repertoire of COVID-19 patients in this study with the 15 CDR3s of antibodies with neutralizing potential collected in two recently published studies. The results revealed that there is a clear preference for sequence homology between clones in this study and antibodies with neutralizing potential. Only the "n3010" antibody from the study of Wu et al has similar sequences to almost all the IGH repertoire in our study, which may reflects the huge diversity of the antibody immune repertoire and the huge differences in the neutralizing activity of antibodies produced by different patients. Patient 4 deserves our attention. According to the results of clone expansion analysis, we can see that the peak time of clone expansion of this patient is much earlier than others. Notably, Patient 4 is the only severe case in this study, and the age is as high as 85 years old. We hypothesized that the earlier peak of clone expansion time point may be assoicated with the severity of clinical symptoms and the age of patients. In general, based on high-throughput immune repertoire sequencing data, this study provides a comprehensive B cell deep immune profiling of COVID-19 patients. Conclusions from this research will help understand the immune response of COVID-19 patients against the SARS-CoV-2 virus and will provide the necessary basic knowledge for the development of vaccines and neutralizing antibodies. A novel coronavirus from patients with pneumonia in China Diversity in the CDR3 region of V(H) is sufficient for most antibody specificities Lineage structure of the human antibody repertoire in response to influenza vaccination Genetic measurement of memory B-cell recall using antibody repertoire sequencing Dynamics of immunoglobulin sequence diversity in HIV-1 infected individuals B-cell repertoire responses to varicella-zoster vaccination in human identical twins Convergent antibody signatures in human dengue Design and standardization of PCR primers and protocols for detection of clonal immunoglobulin and T-cell receptor gene recombinations in suspect lymphoproliferations: Report of the BIOMED-2 concerted action BMH4-CT98-3936 Imonitor: A robust pipeline for TCR and BCR repertoire analysis SOAPnuke: A MapReduce acceleration-supported software for integrated quality control and preprocessing of high-throughput sequencing data IMGT®, the international ImMunoGeneTics information system® Basic local alignment search tool Cd-hit: A fast program for clustering and comparing large sets of protein or nucleotide sequences VDJtools: Unifying Post-analysis of T Cell Receptor Repertoires Clinical Characteristics of Coronavirus Disease 2019 in China Clinical features of patients infected with 2019 novel coronavirus in Wuhan Next-Generation Sequencing of T and B Cell Receptor Repertoires from COVID-19 Patients Showed Signatures Associated with Severity of Disease Lymphopenia an important immunological abnormality in patients with COVID-19: Possible mechanisms High-resolution antibody dynamics of vaccine-induced immune responses High-throughput immunoglobulin repertoire analysis distinguishes between human IgM memory and switched memory B-cell populations The usage of human IGHJ genes follows a particular nonrandom selection: The recombination signal sequence affects the usage of human IGHJ genes Human immune system variation Identification of Human Single-Domain Antibodies against SARS-CoV-2 Potent Neutralizing Antibodies against SARS-CoV-2 Identified by High-Throughput Single-Cell Sequencing of Convalescent Patients' B Cells Network properties derived from deep sequencing of human b-cell receptor repertoires delineate b-cell populations Dynamics of the human antibody repertoire after B cell depletion in systemic sclerosis Immune cell profiling of COVID-19 patients in the recovery stage by single-cell sequencing Deep sequencing of B cell receptor repertoires from COVID-19 patients reveals strong convergent immune signatures The Multifaceted B Cell Response to Influenza Virus Middle east respiratory syndrome coronavirus infection dynamics and antibody responses among clinically diverse patients, Saudi Arabia Broad diversity of neutralizing antibodies isolated from memory B cells in HIV-infected individuals Protective monotherapy against lethal Ebola virus infection by a potently neutralizing antibody. Science (80-. ) Importance of Neutralizing Monoclonal Antibodies Targeting Multiple Antigenic Sites on the Middle East Respiratory Syndrome Coronavirus Spike Glycoprotein To Avoid Neutralization Escape Antibodies and vaccines against Middle East respiratory syndrome coronavirus Antibody 10-1074 suppresses viremia in HIV-1-infected individuals Antibody responses to SARS-CoV-2 in patients with COVID-19 Cross-neutralization of SARS-CoV-2 by a human monoclonal SARS-CoV antibody A noncompeting pair of human neutralizing antibodies block COVID-19 virus binding to its receptor ACE2 Potent human neutralizing antibodies elicited by SARS-CoV-2 infection We sincerely thank the support provided by China National GeneBank. This study was Agency. We thank Dr. Sunil Kumar Sahu at BGI-Shenzhen for his linguistic assistance in this manuscript preparation. The data that support the findings of this study have been deposited into the CNGB Sequence Archive (CNSA: https://db.cngb.org/cnsa/) of CNGBdb with accession number CNP0001106. The authors declare that they have no competing interests.