Differential DNA Methylation and Gene Expression Between ALV-J-Positive and ALV-J-Negative Chickens

Background: Avian leukosis virus subgroup J (ALV-J) is an oncogenic virus that causes serious economic losses in the poultry industry; unfortunately, there is no effective vaccine against ALV-J. DNA methylation plays a crucial role in several biological processes, and an increasing number of diseases have been proven to be related to alterations in DNA methylation. In this study, we screened ALV-J-positive and -negative chickens. Subsequently, we generated and provided the genome-wide gene expression and DNA methylation profiles by MeDIP-seq and RNA-seq of ALV-J-positive and -negative chicken samples; 8,304 differentially methylated regions (DMRs) were identified by MeDIP-seq analysis (p ≤ 0.005) and 515 differentially expressed genes were identified by RNA-seq analysis (p ≤ 0.05). As a result of an integration analysis, we screened six candidate genes to identify ALV-J-negative chickens that possessed differential methylation in the promoter region. Furthermore, TGFB2 played an important role in tumorigenesis and cancer progression, which suggested TGFB2 may be an indicator for identifying ALV-J infections.


INTRODUCTION
As a member of the Alpharetrovirus genus, avian leukosis virus (ALV) causes different pathotypes of neoplastic diseases in chickens (1,2). The spread of ALVs leads to chicken slow growing, production performance degradation and caused serious losses to the poultry industry (3). According to their host range and viral envelope protein, ALVs can be classified as ALV-A, -B, -C, -D, -E, -J, and -K subgroups (4,5). ALV-J was first isolated from meat-type chickens in England in 1988 (4). Since then, ALV-J has been the prevalent subtype of ALVs and has become a serious threat to the world's poultry industry (6,7). To date, there is no effective vaccine against ALV-J. DNA methylation is a major epigenetic mechanism in eukaryotes and plays a crucial role in several biological processes, including the regulation of gene expression, embryonic development, X chromosome inactivation, and the development of various diseases (8)(9)(10)(11)(12). In mammals, DNA methylation occurs mainly at CpG islands and is generally associated with gene repression. However, aberrant methylation has been reported to be associated with various diseases, including neoplastic diseases. A number of studies have indicated that many tumor suppressor genes, such as FHIT, PTEN, and CMTM3, were silenced by promoter hypermethylation in the development of lung cancer, gynecological cancer, and gastric cancer (13)(14)(15). In general, DNA methylation maintains a stable state without environmental stimuli, but many factors can also change the DNA methylation patterns in organisms, including senility, and the development of diet and virus infection (16,17).
The chicken (Gallus gallus) is considered to be an important animal model and, in 2011, the global DNA methylation patterns in chicken genome were analyzed (18). With the development of techniques for sequencing, more and more studies have been designed to identify the genome-wide methylation profiles of chickens, and the results have shown that DNA methylation plays a crucial role in chicken growth development, spermatogenesis, health status, and disease resistant (19)(20)(21)(22). Recently, studies have focused on the interaction between a host response and pathogen infection. Salmonella enterica serovar Enteritidis inoculation has been shown to promote DNA methylation in chicken cecum, which implies that methylation of the HOX gene family may preform important regulatory functions in epigenetic regulation responding to SE inoculation in chickens (23). Marek's disease virus (MDV) can infect chickens causing neoplastic disease and can alter the genome-wide methylation levels of genes in the host (24). The hypomethylation of the CD30 gene occurs in all stages of tumorigenesis, and a high expression of CD30 could be a cue for lymphomas formation after MDV infection (63). ALV-J is also the virus associated with poultry tumors. In our previous study, we examined the aberrantly expressed microRNAs, mRNA, and circRNA in chickens challenged with ALV-J, some of which played important roles in tumorigenesis and the development of post ALV-J infection (25)(26)(27)(28), although genome-wide DNA methylation variation of chickens infected with ALV-J was not fully clarified.
At present, the role of DNA methylation in the pathogenesis of diseases is the subject of intense investigation. In this study, we focused on comparing DNA methylation and gene expression between ALV-J-negative and -positive chickens by performing MeDIP-seq and RNA-seq analyses to identify the differences. Here, we provide basic data and explore the pathogenesis of ALV-J from the perspective of DNA methylation.

Ethics Statement
The animal study protocol was approved by the South China Agricultural University Committee of Animal Experiments (approval ID: SYXK-2019-0136). The experiments were closely followed in accordance with the recommendations of the Guide for the Care and Use of Laboratory Animals of the National Institutes of Health.

Experimental Design
Experimental design is shown in Figure 1. The Chinese local breed chickens were used in this study. All of the 1-dayold chickens were intraperitoneally inoculated with ALV-J NX0101 at a dose of 0.2 mL (10 4.5 TCID 50 /0.1 mL), and inoculated once again at 5 days old. The chickens were raised in negatively pressured biosecurity isolators under quarantine conditions and provided water and feed ad libitum. Whole blood and anti-coagulant blood samples were collected to detect virus or antibodies against ALV-J using reverse transcription PCR (RT-PCR), virus isolation and ELISA assay (29)(30)(31).
Twenty weeks later, chickens tested without ALV-J infection during the period were designated as ALV-J-negative chickens, and the other chickens were identified as ALV-J-positive chickens. The livers of three ALV-J-negative chickens and three ALV-J-positive chickens were selected for MeDIP-Seq and RNA-Seq.

Sample Collection
Three ALV-J-negative chickens and three ALV-J-positive chickens were euthanized to collect the whole livers for further analysis. All the samples were labeled with an ID and transported on dry ice to the laboratory for sample processing and testing. The livers (0.2 g each) were homogenized in phosphate-buffered saline (PBS), and then were detected by Immunofluorescence assay (IFA). Briefly, ALV-J in livers was detected by IFA using the ALV-J-specific monoclonal antibody JE9. Frozen sections of the livers were homogenized, filtered, and then DF-1 cells were inoculated with liver homogenate in a 24-well plate, followed by incubation at 37 • C and 5% CO 2 for 5 days; the cells were incubated with JE9 and detected using FITC-labeled anti-mouse IgG (Biolegend, USA). The fluorescence signals were viewed under a fluorescence microscope (Nikon, Japan). The remaining portion was collected for MeDIP-Seq and RNA-Seq. After this study, the remaining chickens were released to the population for breed conservation.

MeDIP-Seq Library Preparation and Sequencing
Genomic DNA from the samples was isolated using phenolchloroform extraction, precipitated with ethanol, and sonicated to 100-500 bp using a Bioruptor (Diagenode), following the manufacture's protocol (cycle number 6 and cycle conditions (on/off time) 30/30). Sonicated DNA was end repaired, A-tailed, and ligated to adapters by using a NEBNext R Ultra TM DNA Library Prep Kit (NEB). Then, MeDIP-seq was performed with a monoclonal antibody against 5methylcytosine (Active Motif), following the manufacturer's standard protocol. MeDIP DNA libraries were quantified using Quant-iT PicoGreen dsDNA Kits (Life Technologies) and subjected to high-throughput 150 base paired-end sequencing on an Illumina Hiseq sequencer (Cloud-seq Inc., Shanghai, China), according to the manufacturer's recommended protocol.

RNA-Seq Library Preparetion and Sequencing
First, total RNA (1 µg) was used for removing rRNAs using Ribo-Zero rRNA Removal Kits (Illumina, San Diego, CA, USA) following the manufacturer's instructions. Second, RNA libraries were constructed by using rRNA-depleted RNAs with TruSeq Stranded Total RNA Library Prep Kit (Illumina, San Diego, CA, USA) according to the manufacturer's instructions. Libraries were controlled for quality by detecting the library length distribution (Supplementary Figure 1), and quantified using the BioAnalyzer 2100 system (Agilent Technologies, Inc., USA). Lastly, 10 pM libraries were denatured as single-stranded DNA molecules, captured on Illumina flow cells, amplified in situ as clusters, and

Identification of Differentially Methylated Regions
High-quality MeDIP-seq reads were aligned against the Gallus gallus genome (UCSC galGal4) using Bowtie 2 software [v2.2.4, (32)] and only uniquely mapped reads were used for further analysis. Peak calling was performed with MACS software [v1.4.3, (33)]. The peaks in which the midpoint of peaks were located in the region from 2 KB upstream to 2 KB downstream of TSS were defined as promoter peaks; the peaks in which the midpoint of peaks were located in the region from 20 KB upstream to 2 KB upstream of TSS were defined as upstream peaks; the peaks in which the midpoint of peaks were located in the introns were defined as intron peaks; the peaks in which the midpoint of peaks were located in the exons were defined as exon peaks; and the peaks in which the midpoint of peaks were not located in the regions which had been mentioned above were defined as intergenic peaks. Differentially methylated regions (DMRs) were identified by diffReps software (negative binomial test) [v1. 55.4, (34)]. The p-value of the filtering standard was 0.005, and a 2-fold change in the difference of read numbers was used as a criterion for the DMRs.

Identification of Differentially Expressed Genes
RNA-seq high-throughput sequencing and subsequent bioinformatics analysis were done by Cloud-Seq Biotech (Shanghai, China). Briefly, paired-end reads were harvested from the Illumina HiSeq 4000 sequencer, and quality-controlled by Q30. After 3 ′ adaptor trimming and removal of low-quality reads by cutadapt software (v1.9.3), the high-quality clean reads were aligned to the reference genome (UCSC galGal4) with hisat2 software (v2.0.4). Then, guided by the Ensembl gtf gene annotation file, cuffdiff software (35) was used to obtain the gene level FPKM as the expression profiles of mRNA, fold change and p-value were calculated based on FPKM, and differentially expressed mRNA (DEGs) were identified. Genes with a p-value ≤ 0.05 and a log2-transformed value smaller than −1 or greater than 1 were considered to be statistically significant DEGs.

Gene-Ontology (GO) Annotation and the Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway
To further investigate the biological processes and functions associated with differentially expressed genes, we performed GO and KEGG pathway analysis. Genes exhibiting fold change ≥2 and p-value ≤ 0.05 in different samples were analyzed for GO enrichments using clusterProfiler (36,37), and KEGG pathway enrichments using the DAVID functional annotation tool (http:// david.abcc.ncifcrf.gov/). To select the significant GO terms and pathways, a Fisher's exact test followed by the Benjamini-Hochberg (BH) multiple testing correction was performed to calculate the threshold of significance.

MeDIP-Seq Data Validation by Bisulfite Sequencing
To verify differentially methylated genes between ALV-J-negative and -positive chickens, 1 ug of genomic DNA from ALV-Jnegative and -positive chicken samples were treated with sodium bisulfite using an EpiTect Fast LyseAll Bisulfite Kit (QIAGEN, Germany). Primer sequences for the genes selected for validation are documented in Supplementary Table 1.

RNA-Seq Data Validation by Quantitative RT-PCR
A total of 500 ng RNA was used to synthesize cDNA using a PrimeScript RT Reagent Kit (Perfect Real Time, TaKaRa, Osaka, Japan). Quantitative real-time PCR was performed using the 2x SYBR Green qPCR Master Mix (Bimake, USA); qPCR cycling conditions were 95 • C for 5 min, followed by 40 cycles of 95 • C for 10 s, and 60 • C for 30 s. Beta-actin was also amplified and used as a loading control. Relative gene expression was analyzed using the 2 − Ct method. The primer sets used for validation are documented in Supplementary Table 2.

Integrated Analysis of MeDIP and Gene Expression Data
To identify epigenetically regulated genes, the MeDIP-seq and RNA-seq data were integrated. DEGs that retained DMRs in the regulatory regions were selected as candidate genes. Then, genes with upregulated methylated regions (FC ≥ 2, p ≤ 0.005) in the regulatory regions but downregulated expression levels (FC ≥ 2, p ≤ 0.05), or downregulated methylated regions (FC ≥ 2, p ≤ 0.005) in the regulatory regions but upregulated expression levels (FC ≥ 2, p ≤ 0.05) were selected as candidates. All analyses were based on galGal4.

Statistical Analysis
The date from the qPCR shown are mean ± SE from three independent experiments. GraphPad Prism (version 5) was used to process the data. Statistical significance was determined using Student's t-test, and p < 0.05 and p < 0.01 were considered to show significant differences between the groups.

Identification of ALV-J-Negative and -Positive Chickens
To screen the ALV-J-negative and -positive chickens, chickens were monitored for 20 weeks. At 20 weeks, we found that 10 chickens tested negative for ALV-J all the time, and 16 chickens were infected by ALV-J (Supplementary Table 3). Compared to the ALV-J-negative chickens, most of ALV-J-positive chickens showed gradual emaciation. The spleens of the ALV-J-positive chickens were significantly bigger than the ALV-J-negative chickens at 20 weeks (Figure 2A). DF-1 cells were inoculated with liver homogenate from six chickens and the livers of ALV-J-positive chickens were determined to be positive for p27 Frontiers in Veterinary Science | www.frontiersin.org ( Figure 2B). The results showed that we identified the ALV-Jnegative chickens and ALV-J-positive chickens.

Methylomic Profiling of ALV-J-Negative and -Positive Chickens
To screen the methylomic profiling of ALV-J-negative andpositive chickens, the livers of the chickens were sampled, and MeDIP-seq was performed. Following the removal of lowquality data, we obtained a mean of 17,379,530 clean reads from the ALV-J-negative chicken samples and 17,674,372 clean reads from the ALV-J-positive chicken samples. The reads were mapped to the reference genome (UCSC galGal4), and mapping rates of 82.13-83.98% were obtained ( Table 1). The reads were detected across all the mapped chicken chromosomes (Supplementary Figure 2). We analyzed the distribution of peaks among the different genomic components in each sample, including the promoter at the transcription start site, exon, intron, upstream, and the intergenic regions that contained the most peaks. The distribution of methylation peaks in different genomic regions showed similar pattern in those samples (Figure 3).

Characterization of Differential Methylated Regions and Validation of MeDIP-Seq Data Using Bisulfite Sequencing
Differential methylation region (DMR) is considered to be a functional region regulating gene transcriptional level (38).
DiffReps software was used to analyze DMRs in ALV-J-negative and -positive chicken samples. In total, 8,304 DMRs were identified between ALV-J positive and negative chickens. Being that 56.7% of them were hypermethylated and 43.3% were hypomethylated in the ALV-J-positive chickens (Supplementary Table 4). The DMR distribution showed that uniquely mapped reads in intergenic regions had a relatively  higher methylation level than that of other regions, such as promoter, upstream, intron, and exon regions (Figure 4).
To validate the MeDIP-seq data, we mainly focused on the analysis of DMRs located in the promoter region, which may be related to gene expression, therefore, we detected the methylation level of the TGFB2 gene. The bisulfite sequencing PCR (BSP) results confirmed the data accuracy of the MeDIP-seq analysis (Figures 5A,B).
We identified genes that contained DMRs among the methylated genes in ALV-J-negative and -positive chicken samples and GO analysis was carried out using DAVID software. GO terms with p-value ≤ 0.05 were considered to be functionally relevant. The GO analysis showed that 403 terms were enriched, including metal ion binding, DNA binding, and transcription, DNA-templated (Supplementary Table 5, Figure 6). The following four pathways were also identified: steroid biosynthesis, autophagy regulation, beta-alanine metabolism, and ABC transporters (p ≤ 0.05, Supplementary Table 6).

Mapping of RNA-Seq Library Sequencing Data
To identify mRNAs involved in disease resistance, total RNAs from ALV-J-negative and -positive chicken samples were used to construct small RNA libraries. From 22,442,240 to 31,752,754 clean reads were generated from three ALV-J-negative chicken samples and from 19,975,018 to 32,350,232 clean reads were generated from three ALV-J-positive chicken samples ( Table 2). In total, there were 515 differentially expressed mRNAs identified, including 189 upregulated genes and 326 downregulated genes in ALV-J-positive chicken samples as compared with in ALV-J-negative chicken samples (Supplementary Table 7). The heatmap of all samples are shown in Figure 7A. To further understand the function of differentially expressed genes, GO and KEGG pathway analyses were conducted. Most of the GO terms were closely related to the establishment of localization and the cellular development process (Figure 7B and Supplementary Table 8). The KEGG pathway analysis showed that these DEGs were significantly involved in the glycosphingolipid biosynthesis globo series, steroid biosynthesis, and fructose and mannose metabolism pathways ( Figure 7C and Supplementary Table 9). To validate the reliability of the RNA-seq data, qRT-PCR assays were performed. The results showed that qRT-PCR data for six mRNAs were consistent with the observed tendencies using RNA-seq ( Figure 7D).

Integrated Analysis of MeDIP-Seq and RNA-Seq
The integrated analysis between methylation and transcriptome was based on the data from MeDIP-seq and RNA-seq (Supplementary Table 10). The comprehensive distribution patterns of the genes with both differential methylation and expression levels are illustrated in Figure 8A. After merging overlapping DMGs containing DMRs with different gene elements into the unitary DMG, a total of 3,197 DMGs were identified. The genes located in various genomic regions are shown in Figures 8B,C. Furthermore, we observed 72 DEGs that might be regulated by aberrant DNA methylation, including 44 hypermethylation-low expressed genes and 28 hypomethylationhigh expressed genes, which were considered to be potential candidate genes for ALV-J infection. In particular, TMEM104, FAM173A, CNP1, HSD17B7, SUPV3L1, and TGFB2 exhibited differential methylation in the promoter region (Table 3). Furthermore, we also analyzed six of the abovementioned genes with Gallus gallus genome version 6, and the results were similar with those when we used Gallus gallus genome version 4 (Supplementary Table 11).

Bisulfite Sequencing and qPCR Analysis of TGFB2 Expression Level in Chickens
Furthermore, we analyzed the DNA methylation level of TGFB2 in another four randomly selected ALV-J-negative and four ALV-J-positive chickens. As a result, the methylation levels of TGFB2 in ALV-J-negative chickens were higher than that in ALV-J-positive chickens, essentially those results were consistent with MeDIP-seq data ( Figure 9A). In particular,  The DNA methylation and gene expression statuses of the ALV-J-positive chicken samples were utilized as references for those in the ALV-J-negative chicken samples.
there was a methylation site found only in ALV-J-negative chickens. At the same time, we identified whether TGFB2 was commonly downexpressed in ALV-J-negative chickens, the abovementioned chickens which had been analyzed by qRT-PCR. In general, the results were in line with the RNA-seq data ( Figure 9B). Taken together, we found TGFB2 underwent a massive loss of DNA methylation in the promoter region and the expression of TGFB2 was increased during ALV-J infection.

Effect of TGFB2 on ALV-J Replication in DF-1 Cells
In order to understand the reason for TGFB2 with different methylation levels and expression levels, we evaluated the function of TGFB2 in ALV-J replication. We transfected the TGFB2 expression vector into DF-1 cells. The replication of ALV-J was significantly promoted in the pRK5-flag-TGFB2transfected group as compared with the pRK5-flag-transfected group (Figures 10A,B). Meanwhile, to identify whether low expression of TGFB2 would repress replication of ALV-J in DF-1 cells, we transfected DF-1 cells with si-RNA for TGFB2. The replication of ALV-J was significantly repressed in the DF-1 cells transfected with si-RNA for TGFB2, but not si-NC (Figures 10C,D). These results demonstrated the role of TGFB2 in ALV-J replication.

DISCUSSION
Tumor diseases are a serious problem in the poultry industry globally. It is known that ALV-J is associated with several kinds of tumors, such as hemangiomas, erythroblastosis and myelocytomas (39). Vertical transmission of the virus could cause severe immunosuppression in progenies (40), posing significant challenges to the control of this disease. Thus, it is necessary to screen potential molecular markers to identify ALV-J-positive chickens. In this study, we collected liver tissues from ALV-J-positive and -negative chickens and analyzed the DNA methylation because of the tumor emergence in livers of ALV-J-positive chickens in the field. DNA methylation is an essential inheritable modification in most eukaryotic genomes, which plays a crucial role in regulating gene expression. However, this modification is usually altered in tumor cells. The alterations in DNA methylation patterns always affect the susceptibility of cancer. Although several reports regarding the effect of DNA methylation of a gene body on gene expression are available (41,42), the function of DNA methylation of a gene body remains controversial. In contrast, the regulatory effects of DNA methylation of promoter regions on gene expression have been extensively proven in the past years. Estrogen receptor alpha (ERα) which plays an important role in controlling sexual development, is regulated by DNA methylation of promoter regions (43). A change in early growth response 1 (EGR1) expression has been derived from the aberrant methylation of the EGR1 gene promoter region in schizophrenia, and increased ERG1 expression might be connected with the pathophysiology of schizophrenia (44). Concerning glycoprotein VI (GPVI), it has been found that demethylation in the promoter region increased the expression of GPVI, which may be related to the occurrence of coronary heart disease (45). In omental adipocytes, the methylation of CYPl19A1 promoter region has been shown to be negatively associated with relative mRNA expression (46). In patients with type 2 diabetes, DNA methylation of the insulin promoter was increased as compare with non-diabetic donors, and correlated negatively with insulin gene expression in human pancreatic islets (47). All these reports have indicated that the DNA methylation of the promoter regions can regulate gene expression. In this study, we screened over 8,304 DMRs between ALV-J-positive and -negative chickens, and among these, 233 DMRs in the promoter region. More effort is required to analyze the functions of these DMRs in the future.
DNA methylation in the promoter region of genes plays an important role in the intricate host-disease interaction network (48). Abnormal DNA methylation during Marek's disease progression has indicated an interaction between MDV FIGURE 9 | Identification of TGFB2 methylation level and expression in randomly selected chickens by bisulfite sequencing PCR (BSP) and qRT-PCR. (A) The methylation status of TGFB2 in promoter region was investigated among randomly selected chickens by BSP; (B) relative expression of TGFB2 in randomly selected chickens. Data represent means ± SE. **p < 0.01. and the host gene, and has also found that aberrant methylation of the DNMT gene promoter region may also be involved in a virus-induced transformation process (24). The methylation level of ATF5 has been reported to be different in poorly differentiated glioma, well-differentiated glioma, and normal tissues, which suggested that aberrant methylation of ATF5 is connected with the pathophysiology of glioma (49). The methylation levels of SLC26A4 are significantly higher in patients with presbycusis and the methylation at the SLC26A4 promoter can predict the risk of presbycusis (50). Methylation in the p53 promoter region may play a key role in carcinogenesis of epithelial ovarian cancer and has been used as a molecular marker for screening of ovarian cancer (51). SLC5A8 has been shown to be a tumor suppressor in lung tumor and silenced by promoter hypermethylation in the development of lung cancer (52). Several studies have reported that there are some genes with aberrant DNA methylation in promoter regions in endometriosis, such as HOXA10 (53) and ATM (54). Taken together, these reports reveal that DNA methylation is implicated in the development of disease. Until now, there have been few studies regarding the association of DNA methylation with ALV-J infection. We integrated DNA methylation and gene expression data and identified a number of genes which simultaneously changed DNA methylation and gene expression. The integrated analysis showed that the genes expression levels of 72 genes were significantly inversely correlated with DNA methylation level in the ALV-J-negative chickens vs. the ALV-J-positive chickens. These results contributed to the exploration of the potential mechanism of epigenetic studies on the host response to ALV-J infection.
Transforming growth factor β2 (TGFB2), one of the isoforms of TGF-β, has been reported to be associated with various neoplasms (55)(56)(57). TGFB2 can initiate an autocrine loop that promotes its own expression and enables oncogenic activity. In gliomas, the expression levels of TGFB2 are used to evaluate stages of tumor progression (58,59). In ovarian cancer, TGFB2 is overexpressed and plays a key role in ovarian oncogenesis by regulation of an epithelial-tomesenchymal transition (60). TGFB2 has been shown to be changed in different tumor stages, T categories, grades, and patients' survival states, and upregulated in patients with GC as compared with a normal control, and its expression could be affected by cg11976166 (61). In prostate cancer, a quantitative increase in promoter methylation levels of TGFβ2 are associated with PCa progression (62). In this study, ALV-J-negative chickens had a higher level of TGFB2 promoter methylation than ALV-J-positive chickens and the DNA methylation of the TGFB2 promoter had a certain impact on TGFB2 gene expression. In addition, TGFB2 played a significant role in ALV-J replication and was verified in vitro. It is known that ALV-J can cause tumors in chickens; therefore, we suggest that the TGFB2 gene could be a marker gene to identify ALV-J infection.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/, GSE163135 https://www.ncbi.nlm.nih.gov/, GSE143336.

ETHICS STATEMENT
The animal study was reviewed and approved by the South China Agricultural University Committee of Animal Experiments.
Written informed consent was obtained from the owners for the participation of their animals in this study.

AUTHOR CONTRIBUTIONS
YY: conceptualization, methodology, formal analysis, investigation, writing-original draft, writing-review and editing, and visualization. HuiZ: formal analysis, data curation, and visualization. SG: formal analysis, investigation, and visualization. HuanZ: supervision and writing-original draft. XZ: resources. WC: investigation. WL: writing-review and editing and resources. QX: conceptualization, supervision, funding acquisition, and project administration. All authors contributed to the article and approved the submitted version.