H3K4me3 CUT&Tag and Transcriptome Analysis Reveal the Epigenetic Regulatory Landscape in Gill Tissue of Large Yellow Croaker (Larimichthys crocea) Under Low Salinity Stress

H3K4me3 is an important histone modification that could influence DNA replication and RNA translation in response to abiotic stress. Here, RNA-seq analyses were conducted in gill tissues of large yellow croaker to identify the function of H3K4me3 under low salinity stress. Additionally, CUT&Tag analyses were performed to identify the genome-wide dynamic changes in H3K4me3 and explore the mechanisms by which H3K4me3 regulates gene expression. A total of 201 differentially expressed genes (DEGs) were identified between the 5‰ low salinity group (S-group) and 25‰ normal salinity group (C-group), among which 23 DEGs (11 up-regulated H3K4me3 targets and 12 down-regulated targets) were directly regulated by H3K4me3. Our findings thus describe the epigenetic regulatory landscape of H3K4me3 in gill of large yellow croaker during low salinity stress, and provide novel insights into the regulation mechanisms of H3K4me3 mediating the responses of aquatic animals to abiotic stress.


INTRODUCTION
Environmental salinity is a crucial factor that can significantly alter the physiological and endocrine status of marine fish, thereby influencing growth, reproduction, and ultimately survival. Extreme weather events (e.g., heavy precipitation) are becoming increasingly frequent worldwide (Peter et al., 2019). In conjunction with surface runoff, these events can severely affect the salinity of inshore aquaculture ponds and jeopardize the culture of marine organisms. Large yellow croaker is the most productive marine fish artificially cultured in China, and coastal cageculture is the most widely implemented aquaculture model for this species. However, the expansion of large yellow croaker cultivation spaces has become severely restricted due to constant fluctuations in salinity and poor water exchange capacity. Thus, studies on the salinity regulatory mechanisms of large yellow croaker and the selection of varieties capable of adapting to low salinity have recently garnered increasing attention. Most studies on low salinity adaptation of large yellow croaker have mainly focused on ion osmoregulation, energy metabolism, and immune response using omics technology and related molecular biological techniques and methods (Zeng et al., 2017;Lu et al., 2020;Teng et al., 2022). Nevertheless, few studies have assessed the role of epigenetic regulation mechanisms, particularly histone modifications, in marine fish under low salinity stress.
Histone methylation is an important posttranslational modification, which can occur in any components of the nucleosome octamer, but mainly occur on lysine residues, particularly lysine 4 of histone H3 (H3K4), lysine 9 of histone H3 (H3K9), lysine 27 of histone H3 (H3K27), and lysine 36 of histone H3 (H3K36) (Zhang and Reinberg, 2001;Rodriguez and Esteller, 2011). Furthermore, histone methylation levels can also vary (e.g., monomethylation, dimethylation, and trimethylation) and different types of methylated histone can exhibit different biological functions (Bannister et al., 2002;Sims et al., 2003;Esteller, 2007). H3K4 monomethylation (H3K4me1) and H3K4 dimethylation (H3K4me2) sites are more widely distributed than H3K4 trimethylation (H3K3me3) sites, however H3K4me3 plays a critical role in the growth and development of living organisms due to its significantly higher methylation level compared to H3K4me1 and H3K4me2 (Pinon et al., 2017;An et al., 2018). Previous studies have reported that H3K4me3 was highly enriched in the transcription initiation sites (TSS) of genes; and H3K4me3 could both bind with transcription-activated and transcription-inhibited promoters. H3K4me3 binding promoters often have RNA polymerase II binding sites, and therefore the methylation status of H3K4 could significantly influence the transcription initiation of related genes (Mikkelsen et al., 2007;Guenther et al., 2007). Other studies have reported that the methylation status of H3K4 significantly changed under stressful conditions, which in turn induced the expression of stress response-associated genes (Chen and Wu, 2010;Hu et al., 2011;Seibel et al., 2014). H3K4 methylation is a critical driver of the stress response; however, the mechanisms that regulate this phenomenon remain unclear.
Cleavage Under Targets and Tagmentation (CUT&Tag) is a novel method that enables the characterization of the interaction between proteins and chromatin (Hatice et al., 2019). Unlike ChIP-seq, CUT&Tag is highly sensitive due to the high efficiency of tethered Tn5 integration, which also simplifies the library construction and increases the signal-to-noise ratio. This technology thus makes it possible to profile histone modifications and transcription factors on small cell numbers and single cells (Hatice et al., 2020).
Given the quantitative nature of low salinity adaptation, this study comparatively analyzed the dynamic transcriptional and epigenetic changes in gill tissues of large yellow croaker under low salinity stress. By combining CUT&Tag and RNA-seq analyses, this study provides the first comprehensive characterization of transcriptome and epigenome dynamics in marine fish exposed to low salinity stress, thereby enabling a better understanding of the low salinity adaptation mechanisms of large yellow croaker.

Animal Experiments and Sample Collection
Large yellow croaker specimens were collected from aquaculture cages at Ningbo Academy of Oceanology and Fishery in Xiangshan Bay, Zhejiang Province, China. Before the experiment, 150 individuals were temporarily cultured in a 25‰ salinity seawater tank (3-m diameter) for 15 days to acclimate to the experimental environment. The fish was fed daily with commercial mixed feed and the water was exchanged once per day to eliminate residual feed and excrement. Next, 60 healthy similarly sized fish (body weight 76.41 ± 9.55 g, body length 175.01 ± 7.24 mm) were randomly selected and transferred into 6 tanks (1-m diameter; 10 individuals per tank with 500 L 25‰ seawater). After 7 days of acclimation, all of the fish in the 6 tanks were transferred into 6 new tanks respectively which with two salinity treatments [5‰ (S-group) and 25‰ (C-group)] conducted in triplicate. During the experiment, the fish was kept at 22-23°C and DO 8.1-8.7 mg L −1 , fed with commercial mixed feed once per day, and 1/3 of the water was exchanged daily. After 7 days, three individuals from each tank were randomly selected and anesthetized with 50 mg/L MS-222. The gill tissues were then dissected and immediately flash-frozen in liquid nitrogen. The gill samples in each tank were divided into three sections to conduct gene expression, RNA-seq, and CUT&Tag analyses.

Gene Expression Assay
The total RNA of gill tissues was extracted using the TRIzol reagent (Invitrogen, CA, USA), the first strand cDNA was synthesized using the High Capacity cDNA Reverse Transcription Kit (Roache, CA, USA). Quantitative real-time PCR (qRT-PCR) primers were designed using the Primer3 program (Table 1). qRT-PCR was performed using the TB Green Premix Ex Taq II Kit (TaKaRa, Dalian, China). The amplification program consisted of a 95°C pre-denaturation step for 60 s, followed by 40 cycles each of a 95°C denaturation step for 15 s, and a 60°C extension step for 30 s. The expression of target genes was normalized against 18s, which was used as a reference gene.

Transcriptome Sequencing
To characterize the transcriptional dynamics in gill tissues during low salinity exposure, six gill samples from the S-group (5‰) and C-group (25‰) (i.e., triplicate in each group) were collected 7 days after the experimental period. The total RNA of the gill tissues was extracted using the TRIzol reagent (Thermo Fisher Scientifi c, MA, USA). Total RNA concentrations were measured using the Qubit ® RNA Assay Kit coupled with a Qubit ® 2.0 Fluorometer, and its integrity was assessed using an RNA Nano 6000 Assay Kit (Agilent Technologies, CA, USA). The sequencing libraries were constructed using the NEB Next UltraTM RNA Library Prep Kit (New England Biolabs, MA, USA) following the manufacturer's instructions. The libraries were sequenced on an Illumina NovaSeq platform and 150 bp paired-end reads were generated. In order to evaluate the similarity and variability of S-group and C-group, we used DESeq2 software to compare the difference of gene expression among the two salinity group. The clean RNA-seq data were mapped to the reference genome (accession number: GCF_000972845.2). The visual volcano map and DEG heat map was constructed by Ggplot2 software package and pheatmap software package in the R software, respectively. All the differentially expressed genes (DEGs) were mapped to the gene ontology (GO) database and KEGG database using Blast2GO and BLASTx software, respectively. GO and KEGG enrichment analysis was conducted using GOseq R package and KOBAS software, the selection criterion of enriched KEGG pathways was set as P-value < 0.05.

H3K4me3 CUT&Tag Analysis
CUT&Tag analysis was performed using the CUT&Tag 2.0 High-Sensitivity Kit for Illumina (Catalog number: N259, Novoprotein, China) as described by the manufacturers, with some modifications (Hatice et al., 2019). The three gill samples obtained from the S-group (5‰) and C-group (25‰) were then respectively pooled. Briefly, nuclei were extracted from the gill tissues, which were incubated with H3K4me3-specific primary antibodies (AB_2615077) with IgG antibody as a control (normal rabbit IgG, Millipore cat. no. 12-370), followed by the secondary antibody (anti-rabbit IgG antibody, Millipore AP132). A diluted pA-Tn5 adapter complex was then added to the above system, after which DNA was recovered from each of the samples for library construction. Library sequencing was performed on an Illumina NovaSeq platform at Jiayin Biotechnology Ltd. (Shanghai, China).
The FastQC and MultiQC software was used to examine library quality, and the Cutadapt software was used to remove joint sequences and obtain clean sequencing data. Sequence alignment was conducted using the BWA-MEM algorithm of the BWA software, after which sequences were filtered and repetitions were removed using the SAMtools suite. Peak calling and screening were conducted using Macs2 software with a q-value threshold of < 0.05.

Integrated Analysis of CUT&Tag and RNA-Seq Data
To clarify the relationship between H3K4me3 histone modification and related gene expression regulation, we used BETA software (Version 1.0.7) to identify H3K4me3 target genes displaying significant changes that might be attributable to H3K4me3 regulation activity according to a previous study (Wang et al., 2013). Motif analysis was conducted using findMotifsGenome.pl of HOMER software in 200bp up-and down-stream of peaks area; then the predicted motifs were matched to the reported motifs in HOMER and JASPAR database to identify enriched DNA sequences in H3K4me3 binding summits. Further, using the H3K4me3 binding sequence to identify its target genes via comparing with the reference genome of large yellow croaker. Whereafter, the overlapping analysis of the DEGs in RNA-seq and predicted H3K4me3 target genes in CUT&Tag analysis was performed to identify the target DEGs regulated by H3K4me3.

Statistical Analyses
The relative gene expression was calculated using the 2 -△△CT method (Livak and Schmittgen, 2001). Gene expression differences between the C-group and S-group were analyzed via two-way ANOVA followed by Bonferroni's post hoc test.

Expression of Genes Associated With Osmotic Regulation and Histone Methylation
The relative expression of five osmotic regulation and histone methylation related genes (ezh2, prdx3, slc2a10, piezo1, and utx) in the gills of large yellow croaker were calculated after 7 days of low salinity exposure. Our findings indicated that the expression level of these genes were significantly higher in S-group (5‰) compared with C-group (25‰) (Figure 1). Peroxiredoxin (prdx3) is a ubiquitous multifunctional and evolutionarily conserved hyperoxide reductase, which protects cells from reactive oxygen species (ROS)-induced damage (Zhang et al., 2007). Facilitative glucose transporter 10 (slc2a10) is a transmembrane transport protein that mediates the transport of glucose, hexose, and ascorbic acid (Lee et al., 2010). Piezo1 is an evolutionarily conserved electropositive ion transport protein that is involved in cell volume maintenance (Moroni et al., 2018). The changes in the expression levels of prdx3, slc2a10, and piezo1 indicated that the antioxidant level, energy metabolism, and ion transport in the gill of large yellow croaker exhibited adaptive changes during low salinity exposure.

Transcriptional Dynamics in Gill Under Low Salinity Stress
After filtering and removing low quality raw reads, a total of 11.58 million clean reads in the C-group and 11.89 million clean reads in S-group were obtained. For each sample, 6 GB of clean bases were generated. The clean reads of each sample were mapped to the reference genome of large yellow croaker with a 92.45% alignment efficiency, and the ratio of clean reads with unique position in the reference genome exceeded 89.04% (Table 2), with Q30 bases accounting for more than 96% of all samples, confirming that the obtained sequences were highqualified and could thus be used in downstream alignment steps and analyses.
The three replicate samples were grouped together in the cluster analyses, indicating that salinity was the principal stress factor ( Figure 2A). Compared with C-group, a total of 201 DEGs (|log2FC| > 1, p-value < 0.05) were identified in the S-group.
Among these DEGs, 88 were up-regulated and 113 were downregulated ( Figure 2B). Interestingly, more genes were downregulated, indicating that most DEGs were suppressed in gills in response to low salinity stress.
GO terms were categorized as biological process (BP), molecular function (MF), and cellular component (CC). A total of 25,957 unigenes were annotated into 91 sub-categories of the three major GO functional categories, including 36 BP categories, 13 CC categories, and 42 MF categories. The details of gene annotations for the three major categories are illustrated in Figure 2C. The fatty acid metabolic process (GO: 0006631) of the BP category and integral component of plasma membrane (GO: 0005887) of the CC category were significantly enriched. Further, DEGs were significantly enriched into different GO terms with a standard of p-value < 0.02 (Table 3).
Based on DEGs annotation, KEGG analysis provided more detailed information of the biological processes in gill  tissues of large yellow croaker in responses to low salinity stress.

H3K4me3 Modification Changes in Gill Under Low Salinity Stress
In this work, gill samples were collected from large yellow croaker individuals in S-group and C-group after 7 days of exposure. A CUT&Tag analysis was then conducted to evaluate Heatmap of the DEGs between the C-and S-group; red represents up-regulated genes, green represents down-regulated genes, and black represents genes with no significant differences. (B) Volcano plot of the gene expression differences between groups C and S; red represents up-regulated genes, blue represents down-regulated genes, and gray represents genes with no significant differences. (C) GO enrichment of DEGs. The DEGs were annotated into three major functional categories: biological process, cellular component, and molecular function; the x-axis indicates the names of the GO terms and the y-axis indicates the number of DEGs; red represents the total genes between the C and S groups; blue represents the DEGs between the C and S groups. (D) KEGG pathway enrichment analysis. The x-axis indicates the percentage of DEGs in each category; the y-axis indicates the KEGG pathway categories.
the role of H3K4me3 histone modification in low salinity adaptation of large yellow croaker. Qubit analysis was performed after PCR enrichment and purification, and the results indicated that the quality and quantity of the nuclei of the two samples satisfied the requirement of the CUT&Tag reaction. Next generation sequencing (NGS) was then performed, after which 8.53 million and 6.74 million sequences were mapped in the C-and S-groups, respectively ( Table 5). The number of called peaks in the S-group and Cgroup was 2,423 and 6,346, respectively. Further, more than 86% of the reads were mapped to the reference genome for the H3K4me3 histone modification samples. Analysis of the read distribution across genes indicated that the obtained reads count of S-group was significantly higher than that of Cgroup, and the enriched reads, namely the H3K4me3 binding sites, were mainly distributed in the transcriptional start sites (TSSs) (Figures 3A-C). Moreover, most of the reads across peaks in both the S-and C-group were distributed near the peak center, which confirmed the quality and features of the H3K4me3 CUT&Tag data ( Figures 3D-F).
The genomic locations of the peaks were divided into eight categories, including intergenic, intron, exon, 3' UTR, 5' UTR, promoter (<=1kb), promoter (1-2 kb), and promoter (2-3 kb). The H3K4me3 signals were mainly enriched in the intergenic, promoter (<= 1kb), and intron categories in both the S-and C-groups. However, there were significant differences in the percentage of each category between the S-and C-groups ( Figures 4A, B). Computational analysis of H3K4me3 CUT&Tag in large yellow croaker gill tissue under low salinity stress resulted in 1,505 up-regulated peak signals and 4,673 down-regulated peak signals (|log 2 (fold change)| >1, p-value < 0.05). Table 6 shows the location of the lowest 31 peaks ranked by P value. Figures 4C, D show the top five binding motifs of H3K4me3 down-and upregulated targets, respectively. The total significant binding motif number of down-and up-regulated targets was 43 and 39, respectively (enrichment statistics were based on the p-value).
H3K4me3 target genes were predicted via binding and expression target analysis (BETA). A total of 156 predicted H3K4me3 targets were up-regulated and 290 were downregulated, with binding sites within +/-5 kb from TSS. The binding sites of H3K4me3 in large yellow croaker genome were characterized via CUT&Tag. However, this is not sufficient to determine the regulatory landscape of H3K4me3. Thus, an integrative approach using CUT&Tag and RNA-seq data was used to analyze the relationship of H3K4me3 and the target genes expression. The integrated analysis showed that there were 11 up-regulated and 12 down-regulated H3K4me3 target DEGs being identified by overlap the up-and down-regulated DEGs in RNA-seq with the H3K4me3 signal increased and decreased target genes in CUT&Tag, respectively ( Figures 5A, B). GO enrichment analysis indicated that the down-regulated H3K4me3 target genes were associated with signal transduction ( Figure 5C), whereas the up-regulated targets were mainly enriched in the protein kinase C-activating G-protein coupled receptor signaling pathway, intracellular signal transduction, and signal transduction ( Figure 5D). Further, the top 20 KEGG pathways enriched with the target genes were illustrated in Figures 5E, F. The down-regulated targets were mainly enriched in the phosphatidylinositol signaling pathway, inositol phosphate metabolism, T cell receptor signaling pathway, PI3K-Akt signaling pathway, and mTOR signaling pathway. The upregulated targets were mainly enriched in the ether lipid metabolism, glycerophospholipid metabolism, inositol phosphate metabolism, PI3K-Akt signaling, and mTOR signaling pathways.

DISCUSSION
Previous studies have reported that epigenetic modifications could regulate gene expression, genomic rearrangement, and the evolution rate of genomic DNA sequences to acclimate to external environmental pressures (Kim et al., 2005;Asensi et al., 2017). Among the different types of epigenetic modifications, histone modification, DNA methylation, and their interactions have been shown to be major drivers of the acclimatization process (Cedar and Bergman, 2009). Particularly, posttranslational modification of histone plays a vital role in gene transcriptional activation in response to environmental stress, as well as in growth and development. Histone methylation, an important histone modification, is closely associated with gene expression and which is continuously and dynamically regulated by histone methyltransferases (HMTs) and histone demethylases (HDMs). Enhancer of zeste of homolog 2 (ezh2) is a member of polycomb repressive complex 2 (prc2), which possesses histone methyltransferase activity and mediates chromatin transcriptional silencing (Batool et al., 2019). In the present study, ezh2 was significantly activated in the gills of large yellow croaker under low salinity stress, which coincided with changes in histone methylation. Furthermore, utx, a demethylase of H3K27me3, could specifically remove the methylated modification of the lysine residue on H3K27 and promote the transcriptional activation of related genes (Lee et al., 2012). In this study, utx was also significantly up-regulated upon low salinity exposure, indicating that the H3K27me3 histone modification state was also noticeably changed. Therefore, salinity stress could cause a series of histone modification changes in chromatin; however, the functions of these histone modification changes remain unclear.
Studies in other organisms have indicated that H3K4me3 could serve as an indicator of gene activation during the environmental stress adaptation process, and H3K4 methylation could significantly change under environmental pressure (Zong et al., 2013;Singh et al., 2014;Xi et al., 2020). Other studies have indicated that almost 90% of genes in Arabidopsis thaliana carry at least one H3K4 methylated modification under drought stress, and 63% of the genes sites had significant changes of H3K4 trimethylation induced by drought stress (Van et al., 2010). Interestingly, H3K4me3 binding sites are also low standard DNA methylated modification sites, and the interaction of these two epigenetic modifications further enhances chromatin opening and gene transcriptional activation (Du et al., 2015). Very few studies have reported thus far on the functions of H3K4me3 in aquatic animals under environmental stress, and therefore its functional mechanisms remain unclear. The present study provides the first characterization of dynamic transcriptome and epigenome changes in the gills of large yellow croaker under low salinity stress by integrating CUT&Tag and RNA-seq analyses.
The previous studies of low salinity stress were mainly focused on the expression regulation of ion channel and hormone secretion et c. (Larsen et al., 2012;Seale et al., 2014). In this work, GO enrichment analysis indicated that the DEGs in RNA-seq were mainly enriched in metabolism and cellular process and other GO sub-clusters. Thereinto, filamentous actin (F-actin) is the key component of eukaryotic cytoskeleton which plays an important role in cell shape maintenance. Researches indicated that low salinity stress can cause significant morphological changes in gill lamella of large yellow croaker (Teng et al., 2022); the adaptive expression of F-actin can change the structure of fibroid actin microfilament which might be the main cause of morphological change in gill. Phosphopyruvate hydratase is a key enzyme in glucose metabolism which expresses in cytoplasm of almost all organs (Pancholi, 2001). During low salinity adaption, glycolmetabolism related genes' expression in gill could provide adequate energy supply for osmotic regulation, etc. Prefoldin complex is capable of delivering unfolded proteins to cytosolic chaperonin which also plays a vital role in the regulation of chromatin dynamics during transcription elongation (Millań et al., 2013). Besides, low salinity stress has also influenced greatly on erythrocyte homeostasis and homeostasis of number of cells within a gill tissue which is critical in oxygen transport and ion exchange, especially in circulatory system. Previous researches indicated that cell volume recovery under hypotonic stress might rely on calcium signaling pathway; ATP release induced free calcium increase in cell which combined with purinergic receptors in membrane and regulate cell volume recovery (Ollivier et al., 2006). Thereinto, calcium binding proteins is critical in cell volume maintaining, and GTP binding protein has a potential in the regulation of cell calcium concentration (Yuge et al., 2003). Similarly, KEGG pathway analysis indicated that the DGEs enriched in Cell adhesion molecules (CAMs) and Cytokinecytokine receptor interaction pathway were critical in signal transduction and interaction; and the one in tight junction and ferroptosis is potential in maintenance of cell structure and cell apoptosis; besides, the DEGs enriched in PPAR signaling pathway were significantly in expression regulation of related key enzymes in glucose and fatty acid metabolism. Summing up the above, the energy metabolism, signal transduction, and cell function maintenance in gill of large yellow croaker was notably impacted under low salinity stress.
Histone modifications, like H3K4me3, have highly and flexibly response to different environmental stress in the regulation of chromatin remodeling through spatial pattern distribution, and the change of modification condition in specific histone marker genes' promoter region could regulate the expression of environment adaption related genes (Berger, 2007;Zong et al., 2013;Santos et al., 2017). The activation of  drought stress related genes in Arabidopsis thaliana was closely associated with H3K4me3 signal enrichment in gene coding region during drought resistance (Kim et al., 2008). In this work, genome-wide profiling of H3K4me3 data showed that H3K4me3 binding sites were mainly enriched in TSSs, which is consistent with the findings of previous studies (Barski et al., 2007;Fromm and Avramova, 2014), and this distribution pattern is closely related to gene expression and regulation. H3K4me3 could bind to activated and suppressed transcription promoters, thereby up-or down-regulating gene expression (Guenther et al., 2007). Unlike transcript factors (TF) which commonly have unique binding motif (Jairo et al., 2020), H3K4me3 have many binding motifs, and the binding motifs were different from each other between H3K4me3 down-and up-regulated target gene region; this finding gives proof for the flexibility of H3K4me3 regulation to target genes. These dynamic genome-wide H3K4me3 modifications might enable large yellow croaker to adapt to low salinity environments.
The integrative analysis of RNA-seq and H3K4me3 CUT&Tag showed that 11 up-regulated and 12 downregulated H3K4me3 target genes were identified during low salinity exposure of large yellow croaker. According to the GO and KEGG pathway analyses of the H3K4me3 target genes, H3K4me3 modification plays a critical role in signal transduction and energy metabolism, which is consistent with RNA-seq findings. Besides, the H3K4me3 modification condition is closely associated with target genes expression. Unlike H3K46me3, H3K4me3 can both up-and down-regulate target genes' expression. The dynamic changes of H3K4me3 provide flexibility for large yellow croaker in gene's differentially expression under low salinity stress.
In this study, the characterization of transcriptome and epigenome dynamics provided novel insights into the potential mechanisms through which H3K4me3 mediates low salinity adaptation in large yellow croaker. This study thus confirms that histone H3K4me3 modification could participate in the molecular processes that mediate the low salinity adaptation by directly regulating the activation or inhibition of related genes. As to the identified H3K4me3 target DEGs, further research and exploration of the H3K4me3 dynamics and the targets' expression regulation in the low salinity adaptation of large yellow croaker will enable better insights into the underlying H3K4me3 regulating mechanism.

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/, PRJNA798588; https://www.ncbi.nlm.nih.gov/, PRJNA798531.

ETHICS STATEMENT
The animal study was reviewed and approved by Care and Use of laboratory animals of Zhejiang Wanli University. Written informed consent was obtained from the owners for the participation of their animals in this study.

AUTHOR CONTRIBUTIONS
Conceptualization: JY, ZL, and QL. Data curation: JY and TZ. Formal analysis: ML and WS. Funding acquisition: ZL. Methodology: JY. Original manuscript writing: JY. Manuscript review and editing: ZL and QL. All authors contributed to the article and approved the submitted version.