Chromatin accessibility and differentially expressed genes profiling in large yellow croaker (Larimichthys crocea) head kidney cells following iridovirus infection

Abstract

Introduction:

The large yellow croaker iridovirus (LYCIV) poses a significant threat to the aquaculture industry of Larimichthys crocea. Understanding the host defense response to LYCIV infection is crucial for developing effective strategies to mitigate its impact.

Methods:

In this study, an epigenetic approach was employed to investigate dynamic changes in chromatin accessibility using the assay for transposase-accessible chromatin with high-throughput sequencing (ATAC-seq). Additionally, RNA sequencing (RNA-seq) was used to analyze the expression pattern of immune response genes upon LYCIV infection.

Results:

Substantial alterations in chromatin accessibility were observed, particularly in the regulatory regions of key immune-related genes. Significant changes in the expression of AP-1 transcription factors, including the Batf gene, were noted. CUT&Tag results revealed that AP-1 was significantly enriched in the open chromatin regions of cytokine genes, with Batf potentially regulating the cytokine genes LIF and CLCF1.

Discussion:

These findings suggest that AP-1 may play a crucial role in the defense response against viral infection by modulating inflammatory cytokines and contributing to cellular inflammatory responses. This study provides a comprehensive analysis of the epigenomic landscape and gene expression regulation during iridovirus infection in L. crocea, offering valuable insights for breeding programs aimed at combating iridovirus infections.

1 Introduction

The large yellow croaker (Larimichthys crocea) with its distinctive golden yellow body color is of great economic importance in marine fisheries in China. According to the 2023 China Fisheries Statistical Yearbook, nationwide aquaculture production of large yellow croaker reached 257,700 tons in 2022, making it the leading marine aquaculture species (1). The iridovirus disease is a prevalent viral infection in large yellow croaker. Iridovirus is an enveloped dsDNA virus of the Iridoviridae family, with a genome size of about 112 kb (2, 3). According to the International Committee on Taxonomy of Viruses (https://ictv.global/taxonomy), iridoviruses are currently classified into seven genera, namely Iridovirus, Megalocytivirus, Ranavirus, Lymphocystivirus, Chloriridovirus, Decapodiridovirus, and Daphniairidovirus (4). Among these, Megalocytivirus represents one of the most significant pathogens responsible for high mortality in many finfish species (5). The genus Megalocytivirus comprises two species, Megalocytivirus lates1 and Megalocytivirus pagrus1. The large yellow croaker iridovirus (LYCIV) belongs to the Megalocytivirus pagrus1 and is closely related to red sea bream iridovirus (RSIV), infectious spleen and kidney necrosis virus (ISKNV), turbot reddish body iridovirus (TRBIV) and rock bream iridovirus (RBIV) (6), all of which are part of the Megalocytivirus genus (4). The infection normally manifests as body darkening, reduced appetite, and severe tissue damage in organs including the liver, spleen, and kidneys, with mortality rates reaching up to 75% (7). Understanding the defense response of large yellow croaker against iridovirus attack, especially the induced-gene expression profile, is crucial for advancing knowledge in virology and immunology (8).

An interesting topic of epigenetic regulation of gene expression is the impact of chromatin structure and organization. The structural state of chromatin directly influences the interactions between transcription factors and the genome, facilitating the access of transcription machinery (9). In other words, the chromatin accessibility is a pivotal switch for the potential onset of regulated gene expression, with the open chromatin region favoring the initiation of regulatory activities (10). The accessibility of these open chromatin regions is not static but dynamic (11). This dynamic nature allows genes within the genome to appropriately respond to external stimuli, facilitating processes such as cell differentiation and defense responses (12).

Assay for Transposase Accessible Chromatin with high-throughput sequencing (ATAC-seq) is a powerful tool to probe chromatin accessibility at whole genome level (13). It has been applied to investigate changes in host chromatin accessibility following pathogen infection and their effects on immune responses (14). The integration of ATAC-seq and RNA-seq identified candidate transcription factors (TFs) involved in viral responses in rice infected with RSV (15). The comparison of chromatin accessibility in phagocytic versus non-phagocytic cells in the head kidney of grass carp revealed that increased phagocytic capacity was correlated with reduced chromatin accessibility at promoter regions, which affected the accessibility of the transcription factors Irf4 and Irf8 crucial for B cell development (16–18). Moreover, analysis of macrophage chromatin state in obesity model mice versus control mice identified ETV5 as a critical transcription factor (19). These applications indicate a great potential of ATAC-seq as a useful technique to uncover the host immune regulatory network by epigenetic modulation.

Viral infections can significantly increase the chromatin accessibility of host antiviral regulatory elements, thereby enhancing the expression of antiviral genes. The scaffold protein SAFA (hnRNPU) can regulate the opening and closing of antiviral regulatory elements to regulate the expression of antiviral immune genes upon recognition of viral double-stranded RNA presence (20). Disruption of SAFA expression in cells notably reduces the chromatin accessibility of virus-inducible genes, while the accessibility of non-virus-responsive genes remains largely unaffected (21). On the other hand, some DNA virus infections can change the location and organization of host chromatin leading to an extensive remodeling to facilitate the viral replication (11). Investigation of host chromatin structural changes in accessibility across the genome in response to viral infection is important for a better understanding of host-virus interaction and host immune response. Currently, there is limited research on chromatin accessibility changes in fish following viral infection, highlighting the necessity of further exploration in this field.

In this study, the ATAC-seq technology was employed to study the dynamic changes in chromatin accessibility in large yellow croakers following iridovirus infection, and the coupled RNA-seq investigation was performed to explore how these changes influence the expression of host immune genes. Analysis of ATAC-seq and RNA-seq data on large yellow croaker cells at various time points of iridovirus infection revealed that the viral infection significantly altered chromatin openness, particularly in the regulatory regions of key immune-related genes. These findings offer new insights into the intricate antiviral activities and contribute to the development of novel antiviral strategies.

2 Materials and methods

2.1 Experimental fish collection and virus identification

About 100 juvenile large yellow croakers with typical symptoms of iridovirus infection were collected from an aquaculture center (Ningde City, Fujian Province), for the preparation of a virus infection solution. Diseased fish were identified by observing the characteristics of diseased fish (body darkening, reduced appetite, and severe tissue damage in organs including the liver, spleen, and kidney) and using specific PCR primers. These symptomatic fish had an average body length of 9.57 ± 2.43 cm, a body height of 2.52 ± 0.48 cm, and a body weight of 11.25 ± 1.75 g. The fish harvested were immediately frozen in liquid nitrogen and stored at -80°C to maintain sample integrity. The sample collection and experimental procedures were approved by the Jimei University Animal Care and Use Committee and strictly adhered to the “Regulations for the Administration of Affairs Concerning Experimental Animals” issued by the State Council of China. Virus suspension preparation was performed as described in previous studies as following (3). To prepare the virus infection solution, the liver, spleen, and kidneys of the diseased fish after thawing were harvested in sterile 2 mL microcentrifuge tubes and homogenized on ice using a tissue homogenizer (60 Hz, 120 s). After dilution in serum-free DMEM (1:10 w/v) supplemented with double antibiotics (penicillin-streptomycin) and amphotericin B, the homogenized samples were incubated in a water bath at 15°C for 2-4 h or at 6°C for 6-24 h. The samples were then centrifuged at 4°C at 5,000 rpm for 5 min. The obtained supernatant was filtered through a 0.22 µm filter and the filtrate containing virus was stored at -80°C.

Viral DNA was extracted using a specialized viral DNA extraction kit (ER201-01, TransGen). PCR amplification was performed on the extracted DNA using LYCIV-specific primer pairs LYCIV-F1 and LYCIV-R1 (Table 1) targeting the major capsid protein (MCP) gene (GenBank OK149102.1), and LYCIV-F2 and LYCIV-R2 (Table 1) targeting the ATPase gene (GenBank KY774434.1), respectively. The amplicon was examined by agarose gel electrophoresis.

2.2 Viral load quantification and cell culture

A 197-bp PCR amplicon by the primer pair of LYCIV-F3/R3 (Table 1) targeting the iridovirus ATPase gene (Supplementary Figure S1) was constructed into the pEASY-T1 cloning vector. The recombinant plasmid was transformed into Trans1-T1 competent cells for replication. After isolation and quantification, a serial 10-fold dilution of the recombinant plasmid was prepared. A viral load standard curve (R²= 0.9981) was established using the serial dilutions as the standards by real time quantitative PCR (RT-qPCR) (Supplementary Figure S2), which was used as a reference for iridovirus viral load determination based on the cycle threshold.

The large yellow croaker head kidney (LYCK) cell line (22) was cultured in Leibovitz’s L-15 medium supplemented with 10% fetal bovine serum at 28°C. Cells in optimal growth conditions were digested with trypsin and seeded into 6-well plates. Once the cells reached 80% confluence, they were infected with iridovirus using an infection solution with a concentration of 5×10⁴ copies/µL and a multiplicity of infection (MOI) of 100. The cultures were then maintained at 28°C and sampled at 0, 12, 24, and 36 hours post infection (hpi), respectively, with three replicates at each time point. After centrifugation to remove trypsin, the cells were resuspended in 10 μL of PBS for subsequent library construction for ATAC-seq or RNA-seq.

2.3 ATAC-seq library preparation and data analysis

After determining the total and viable cell count in the cell resuspension using a Luna II cell counter, approximately 50,000 viable cells were harvested by centrifugation at 800 g for 5 min. The cells were resuspended in 50 μL of precooled lysis buffer containing 10 mM Tris-HCl (pH 7.4), 10 mM NaCl, 3 mM MgCl₂, and 0.1% IGEPAL CA-630, and lysed on ice for 10 min to release the nuclei. The nuclei were collected by centrifugation at 500 g for 5 min in 4°C and subjected to a transposition reaction and library construction using the TD501 ATAC-seq Kit (Vazyme, China). After the transposition reaction, DNA was captured and purified using pre-equilibrated magnetic beads, followed by two washes with 80% ethanol. Library fragment size distribution was ensured to match the expected mononucleosome or dinucleosome model using Qseq-100 (Bioptic, China). The total concentration of the library was estimated using Qubit 3.0 fluorometer (Life Technologies, US), and high-throughput sequencing was performed on Illumina NovaSeq platform with a 150 bp paired-end mode (23).

After data trimming and quality control using Trim-galore (v2.10) (the parameters are set to -q 25 –phred33 –length 35 -e 0.1 –stringency 4 –paired) (24), and FastQC (v0.11.9) (25), the clean reads were aligned to the reference genome of L. crocea (https://www.ncbi.nlm.nih.gov/assembly/GCF_000972845.2#) using Bowtie2 v7.5.0 (26), the parameters are set to -p 5 –very-sensitive -X 2,000. The duplicate reads, mitochondrial reads, and low-quality alignment reads were filtered out using Sambamba v0.7.1 to obtain high-quality alignments (27) using default parameters. The enrichment of transcriptional initiation sites (TSS) was identified using DeepTools v3.4.2 (28). Peaks were identified using Macs2 v2.2.7.1 (29) (-g 7.5e08 –nomodel –shift -100 –extsize 200). The consensus peaks from three replicates of each stage were determined using Bedtools v2.29.2 (30). Annotation of genomic locations of the ATAC-seq peaks was achieved using ChIPseeker v1.24.0 package (31). Significantly enriched motifs in peaks were identified using the tool of Homer v4.10 (32), with motifs of 8, 10, and 12 in length.

To analyze open chromatin regions, the Diffbind package (33) was used to calculate read counts for each peak, adjust peak intervals between samples, centralize and prune peaks based on read overlaps, and evaluate differences based on edgeR with a significance threshold of FDR < 0.05 (34). The 5 kb in front of TSS was defined as the potential promoter region, and the position of peaks near each gene was distinguished, such as intron region, distal intergenic region, and so on. An online bioinformatics tool jvenn (https://www.bioinformatics.com.cn/static/others/jvenn/example.html) was used to obtain significant enrichment motifs intersections of peaks at each infection point, and those common motifs were identified.

2.4 RNA-seq library construction and data analysis

Total RNA was extracted using TRIzol reagent (Invitrogen, Carlsbad, CA, United States), and its quality was assessed by agarose gel electrophoresis and further analyzed using a Nanodrop spectrophotometer and a Qubit fluorometer and Agilent 2100 bioanalyzer. The library was created using TIANSeq Stranded RNA-Seq Kit (Illumina) and subjected to the quality control analysis using Qubit fluorometer and Agilent 2100 bioanalyzer. Once the inserted fragment size met our expectations, the library was quantified using qPCR. The libraries were then pooled and sequenced on the Illumina Novaseq 6000 sequencing platform.

The RNA-seq data was cleaned using STAR v2.7.6a (35) to align clean reads to the reference genome (https://www.ncbi.nlm.nih.gov/assembly/GCF_000972845.2#), and the transcripts were assembled using StringTie v2.1.2 (36). Transcript abundance was calculated using featureCounts v1.0 (37), and differentially expressed genes (DEGs) were identified using DESeq2 v1.28.1 (38). The criteria for defining differentially were |log2(fold change)| > 1 and adjusted p-value < 0.05. The Gene Ontology (GO) enrichment analysis of the DEGs within each group was performed using clusterProfiler v3.16.1 (39). The AnnotationHub Server (v2.20.2) (40) was used to obtain GO annotation information for the large yellow croaker. The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment was analyzed with DOSE v3.14.0 (41), with a p-value < 0.05 for significant enrichment.

SYBR Green-based RT-qPCR was performed to analyze the expression of interleuin-1beta, irf8, il13ra1-like, and ccl19 (genes involved in the cytokine-cytokine receptor interaction pathway) using the 2^−ΔΔCT method. RT-qPCR was also performed to analyze the expression of fgf21, g-lysozyme, pim1, rpa1, sbspon and cdk6 (genes involved in processes such as cellular stress, DNA damage response, and immune regulation) using the 2^−ΔΔCT method. β-actin was chosen as the reference gene for normalization (see primers in Table 1).

2.5 Batf gene expression analysis

An overexpression vector was designed based on the CDS sequence of the batf gene (NCBI: XM_010735250.3), and primers were designed and synthesized by Sangon Biotech (Shanghai, China). The single-stranded primers were annealed into double-stranded Oligo sequences and connected to the double-enzyme linearized overexpression vector (pBHGlox(delta)E13Cre). The correct transformants were verified by sequencing, and high-purity plasmid extraction was performed. Thus, a plasmid containing the Ad-shuttlescr-batf overexpression vector was obtained. The adenovirus shuttle plasmid carrying the exogenous gene and the auxiliary packaging plasmid carrying most of the adenovirus genome were co-transfected into HEK293 cells, and the recombination was achieved using the Cre/loxP recombinase system to produce recombinant adenovirus (38, 39). The HEK293T and LYCK cells were transfected with a negative control adenovirus vector and those transfected with the Batf-overexpression adenoviral vector. RT-qPCR was performed to analyze the induced expression of batf gene with β-actin as the internal control (Table 1). Western blot was performed following the standard procedures using anti-Flag as the primary antibody.

2.6 Adenovirus infection and CUT&Tag library construction

LYCK cells were detached and diluted to 3×10⁵ cells/mL. Subsequently, 1.5×10⁵ cells (500 μL/well) were seeded into each well of a 24-well plate and incubated overnight at 28°C. Adenovirus was thawed on ice and adjusted to viral titer of 1.07×10¹¹ PFU/mL at the time of infection. Based on preliminary experiments, an MOI (multiplicity of infection) of 100 was deemed suitable. After 4 h of infection, 500 μL of complete culture medium was added. Following an additional 6-8 h of incubation, the virus-containing culture medium was removed, replaced with fresh complete medium, and cells were further cultured at 28°C.

The cells were detached using trypsin digestion. Once completely dissociated from the culture surface, the digestion was halted using complete culture medium. The cells were then mixed, centrifuged at 800 g for 5 min. After removing the supernatant, 1 mL of PBS was added, and a 1:1 mixture of 10 μL of cell suspension and Trypan blue was prepared, and the cells were counted using a cell counter with three biological replicates. The cultured cells were then washed and bound to ConA beads. The cell-ConA bead complex was then co-incubated with the primary antibody, followed by incubation with the secondary antibody. Hyperactive pA/G-Transposon was added for TTBL tagmentation. The DNA was then extracted, and the library was amplified and purified. The distribution and concentration of library fragments were detected (39, 40).

The parts containing Adapter, N, or low-quality bases in the raw sequencing data were filtered out by Cutadapt v4.8.137 (42) to obtain high-quality sequencing data (clean reads) for follow-up alignment. Bowtie2 v7.5.0125 (26) was used to align the high-quality sequencing data with the reference genome (https://www.ncbi.nlm.nih.gov/assembly/GCF_000972845.2#), and then the results were evaluated for better quality. After removing reads from the mitochondrial sequence, the CUT&Tag signal (peaks) was scanned across the entire gene range using peak-calling software Macs2 v2.2.7.1 (29). The position information of peaks in the genome and the Peaks length information were obtained. Bedtools v2.29.2 (30) was used to merge two biologically consistent Peaks.

2.7 Data availability

The datasets [ATAC-seq, RNA-seq, CUT&Tag] generated in this study can be found in the [PRJNA820167], and [PRJNA1116991] of Sequence Read Archive (SRA).

3 Results

3.1 Identification of iridovirus and cytopathic effects

The presence of large yellow croaker iridovirus (LYCIV) in various organs (liver, spleen, mid-kidney, gills) of the symptomatic large yellow croaker was confirmed by PCR. A specific band was detected in the diseased fish but was absent in the non-template control (Figure 1A). In addition, the LYCIV was specifically detected in the large yellow croaker cell line and the culture media post infection (Figure 1B), which indicated a successful colonization of LYCIV in the LYCK cells.

Figure 1

Phenotypically, the LYCK cells post infection exhibited noticeable vacuolization and other cytopathic effects (CPE), which was intensified over the course of the infection, including increased cellular vacuolization and contraction (Figure 1C).

3.2 Peak calling and genomic distribution of ATAC-seq peaks

The ATAC-seq library at each time point of infection was constructed and rigorous quality control was performed. After analysis of the ATAC-seq data, the clean reads, the alignment rate to the reference genome, and the enrichment at transcriptional starting sites (TSS) were listed (Table 2). The average size of the ATAC-seq clean data was 7.98 Gb, and the average alignment rate to the reference genome (GCF_000972845.2) was above 90% for each group, suggesting a good quality of the sample preparation. The length distribution of insert fragments showed three distinct nucleosome models, with peaks at approximately 190 bp, 350 bp, and 520 bp, respectively (Figure 2A). TSS enrichment analysis revealed significant enrichment of signals near transcription start sites (Figure 2B). Most of the ATAC-seq peaks were located in gene promoter regions, introns, and distal intergenic regions (Figure 2C).

Figure 2

3.3 Dynamic changes in chromatin accessibility and motif analysis

Systematical analysis of chromatin accessibility changes in LYCK cells in response to LYCIV infection showed that 12,201 regions exhibited decreased chromatin accessibility, while 8,884 regions exhibited increased accessibility at 12 hpi compared to the control (P12_VS_P0). Additionally, 16,167 regions had decreased and 16,228 regions had increased accessibility at 24 hpi compared to those at 12 hpi (P24_VS_P12). Furthermore, 325 regions showed decreased and 181 regions showed increased accessibility at 36 hpi compared to those at 24 hpi (P36_VS_P24) (Table 3). The dynamic changes of chromatin accessibility suggested an active chromatin remodeling in LYCK cells in response to the viral attack. The percentage of upregulated ATAC-seq peaks in the promoter regions increased consistently, while the percentage in distal intergenic regions gradually decreased (Figure 2D). The top 20 enriched motifs in regions with chromatin accessibility changes across the three comparisons (P12_VS_P0, P24_VS_P12, P36_VS_P24) were analyzed and showed in the Venn diagram (Figure 2E), and some motifs common to all three time points were identified such as AP-1 transcription factor motifs. The AP-1 family involved includes Batf, JunB, Atf3 and other transcription factors, which shared a core binding motif of TGACTCA (Figure 2F).

3.4 RNA-seq analysis and pattern of DEGs

After RNA-seq library preparation, the quality of the library was assayed using Agilent bioanalyzer. Based on the strip sizes in the simulated gel chart, 28S, 18S, 5S, and Lower Marker bands were identified (Figure 3A). All 12 libraries passed quality control and met the sequencing requirements. RNA-seq of the 12 libraries produced a total of 85.22 Gb of data, with an average of approximately 7.10 Gb per library. After filtration, 81.85 Gb of data remained, with an average of about 6.82 Gb per library. The Q20 ratio ranged from 95.38% to 97.64%, the Q30 ratio from 90.00% to 93.86%, and the GC content from 44.46% to 47.77% (Table 4).

Figure 3

In the P12_VS_P0 group, 1,089 genes were significantly downregulated, and 1,002 genes were significantly upregulated (Figure 3B). In the P24_VS_P12 group, 84 genes were significantly downregulated, and 140 genes were significantly upregulated (Figure 3C). In the P36_VS_P24 group, 1,323 genes were significantly downregulated, and 2,094 genes were significantly upregulated (Figure 3D). The number of differentially expressed genes at each time point was shown in the chart (Figure 3E).

3.5 GO and KEGG analyses of DEGs

Gene ontology (GO) analysis revealed that the upregulated gene set in P12_VS_P0 (P12_VS_P0_UP) was significantly enriched in pathways related to DNA damage response (Figure 4A). Concurrently, the downregulated gene set in P24_VS_P12 (P24_VS_P12_DOWN) showed significant enrichment in pathways associated with immune system regulation (Figure 4B). The upregulated gene set in P36_VS_P24 (P36_VS_P24_UP) was notably enriched in processes related to negative transcriptional regulation (Figure 4C). In contrast, the downregulated gene set in P36_VS_P24 (P36_VS_P24_DOWN) was significantly enriched in pathways related to DNA damage response, stress response processes, and pressure response regulation (Figure 4D).

Figure 4

KEGG analysis showed the gene sets in P12_VS_P0_UP and P24_VS_P12_DOWN was significantly enriched in the cytokine-cytokine receptor interaction pathway. Similarly, the gene set in P12_VS_P0_UP and the gene set in P36_VS_P24_DOWN were significantly enriched in the lysosome pathway (Figure 4E).

To validate the RNA-seq data, a set of genes were selected to examine their expression pattern by RT-qPCR, including the genes involved in cellular stress, DNA damage response, and immune regulation. The analysis of transcript abundance represented by TPM (transcripts per million) in RNA-seq and the relative expression level represented by normalized RT-qPCR showed a good correlation, indicating good reliability of the RNA-seq data (Figures 5A–J).

Figure 5

3.6 Synergistic changes in chromatin accessibility and gene expression

To further elucidate the dynamic changes in chromatin accessibility and gene expression during LYCIV infection in LYCK cells, the differential chromatin accessibility peaks from ATAC-seq and the DEGs from RNA-seq were compared with different combinations. The results indicated that the number of upregulated peaks was greater than that of downregulated peaks near upregulated genes, while the number of downregulated peaks near downregulated genes was greater than that of upregulated peaks, indicating a positive correlation between the chromatin accessibility and gene expression in a good amount of genome regions (Figure 6A). For example, the promoter region of the upregulated cytokine gene il1 exhibited increased chromatin accessibility after 12 h of infection, as evidenced by a substantial increase in ATAC-seq peaks. Concurrently, the gene expression significantly increased, as indicated by the higher stacking height of RNA-seq reads on the transcript exon (Figure 6B).

Figure 6

Some key transcription factors that showed significant upregulation of gene expression and significant enrichment in chromatin accessible regions following infection were of great interest. Among the four common transcription factors (BATF: AP-1, MAFK, BACH2, and NFATC3), BATF, a member of the AP-1 transcription factor family, was significantly upregulated with a p-value < 0.01 in the comparison between P12 and P0 of RNA-seq data (Figure 6C). Simultaneously, there was significant enrichment of the BATF motif in the upregulated chromatin accessibility regions, with a p-value < 0.01 in the comparison between P12 and P0 (Figure 6D), further substantiating the relationship between chromatin state and regulatory gene expression.

From the non-coding regions of differentially expressed cytokine genes (promoter, intron, and intergenic regions), significantly enriched transcription factor binding motifs were extracted from the ATAC-seq peaks. The results showed that the TRE (TGAC/GTCA) regulatory element of AP-1 was significantly enriched in the open chromatin regions of these non-coding regions of cytokine genes (Figure 6E).

3.7 Target verification of CUT&Tag

To verify whether the key transcription factor AP-1 regulates the expression of cytokine genes, Batf gene was cloned into an adenovirus vector with GFP expression cassette and transfected into LYCK cells. Fluorescence microscopy showed strong expression of Batf-GFP (Figure 7A), suggesting successful transformation and expression of Batf gene. In addition, the expression of Batf gene was examined by RT-qPCR and western blot. The RT-qPCR indicated abundant transcripts of Batf gene (Figure 7B), and the Flag-tagged Batf protein was detected specifically in the transfected cells with the Batf expression cassette (Figure 7C).

Figure 7

The antibody-guided CUT&TAG library preparation showed that the H3K4me3 library exhibited a regular pattern of three nucleosome models, indicating correct library preparation. In contrast, the experimental Batf and IgG libraries displayed irregular patterns, with the IgG library overall showing lower concentrations (Figure 7D). The offline data from the CUT&Tag sequencing were filtered and aligned to the reference genome of L. crocea, achieving an alignment rate of over 90%, with an average of 91.87%, which was sufficient for subsequent analysis (Table 5).

The CUT&Tag assay results demonstrated that Batf bound to the open chromatin regions of the cytokine gene LIF. Peaks of H3K4me3 and open chromatin were observed in the promoter region of the LIF gene, along with Batf-targeted binding peaks, but no signal from IgG (Figure 7E). Similarly, Batf was bound to the open chromatin region of the cytokine gene CLCF1 (Figure 7F). In the promoter region of the CLCF1 gene, peaks of H3K4me3 and open chromatin were present, along with Batf-targeted binding peaks, and no signal from IgG. These findings were consistent with the results from previous ATAC-seq and RNA-seq analyses, further confirming the targeted regulation of cytokine genes by AP-1.

4 Discussion

The results of ATAC-seq revealed significant changes in chromatin accessibility within open chromatin regions during the progression of iridovirus infection. These changes likely reflect how host cells regulate gene expression at the chromatin level to adapt to iridovirus infection. Iridovirus infection induced significant alterations in chromatin architecture, leading to the reorganization of host gene expression patterns. This dynamic change in chromatin accessibility is comparable to the changes observed in host chromatin accessibility caused by Bombyx mori nucleopolyhedrovirus (BmNPV) infection in previous studies (43). We observed that iridovirus infection resulted in a notable increase in chromatin accessibility at several immune-related loci, particularly in non-coding regions regulating the cytokine-cytokine receptor interaction pathway. This is consistent with previous reports indicating that viral infections can induce changes in host chromatin structure to enhance the expression of immune genes (20). For instance, the scaffolding protein SAFA (hnRNPU) has been shown to recognize viral double-stranded RNA and selectively modulate the accessibility of antiviral regulatory elements, thereby boosting the expression of immune genes (20).

In our study, the significantly enriched GO and KEGG pathways of differentially expressed genes at various time points suggested a response pattern of LYCK cells to iridovirus infection. This view is supported by KEGG analysis and aligns with previously reported results of host defense failure caused by rock sea bream iridovirus (RBIV) (44).

We found that the interaction between cytokines and cytokine receptors is the most significantly enriched pathway in the KEGG analysis. However, one factor we cannot exclude is that regulatory elements responsible for upregulating chromatin accessibility might also be associated with the expression of non-cytokine genes. Additionally, it is possible that the upregulation of cytokine gene expression could be related to regions where chromatin accessibility remains unchanged nearby. To address these concerns, we conducted motif analysis on the open chromatin regions of all cytokine genes to identify upstream transcription factors regulating cytokine gene expression. This analysis revealed significant enrichment of AP-1 motifs in the open chromatin regions of non-coding regions associated with cytokine genes, suggesting a potential targeted regulatory role of the transcription factor AP-1 on cytokines.

ATAC-seq is highly sensitive in detecting potential transcription factor binding sites, offering valuable insights into the regulatory mechanisms activated in response to iridovirus infection (45). Our study identified AP-1 as a key transcription factor, with its binding sites showing increased chromatin accessibility following infection. Activator Protein-1 (AP-1) is crucial for various biological processes. It functions as a heterodimer composed of members from different DNA-binding protein families. AP-1 binds to TPA response elements (TRE) and cAMP response elements (CRE; TGAC/GTCA) (46). Its diverse functions are part of a complex and dynamic network of signaling pathways, which depend on the composition of its subunits and their interactions with other nuclear factors (47). In this study, we observed that the batf gene, a member of the AP-1 family (basic leucine zipper ATF-like transcription factor), was significantly upregulated 12 hpi and was notably enriched in regions of increased chromatin accessibility.

Nevertheless, further validation of AP-1’s targeted regulation is required. Conducting AP-1 ChIP-Seq experiments would provide more definitive evidence; however, performing ChIP-Seq in non-model species presents challenges such as high background noise, poor reproducibility, lack of highly specific antibodies, and the complexity of library construction (48). To overcome these limitations, our study employed adenoviral infection of LYCK cells to introduce an overexpression of Batf and take advantage of the tagged protein for CUT&Tag assay. In the western blot experiment involving the overexpression of Batf, a band near 55 kDa was observed, which may represent the dimer form of the Batf transcription factor. Previous studies have demonstrated that the leucine zipper domain of Batf facilitates dimerization with members of the Jun protein family (49). Additionally, some bands near 25 kDa were also observed in the western blot assay for Batf overexpression. These bands may be related to the phosphorylation of JNK kinase, which is associated with the activation of AP-1 transcriptional binding activity (50). AP-1 members frequently undergo post-translational modifications, such as phosphorylation by JNK kinase. The addition of phosphate groups increases the negative charge of the phosphorylated protein, leading to slower migration through the gel and an apparent increase in molecular weight, which was observed in c-Fos, another member of the AP-1 family (51).

In CUT&Tag experiment, the transposase cuts at the regions where the primary and secondary antibodies are bound. H3K4me3 (trimethylation of lysine 4 on histone H3) is typically located in open chromatin regions, so the transposase-antibody complex that binds to H3K4me3 will cut nearby nucleosomes, resulting in a characteristic nucleosome pattern in the library length distribution showing a regular pattern of three nucleosome models in a well-prepared library. On the other hand, IgG binds non-specifically to genomic regions and can be used to exclude non-specific binding regions in the experiment. Results from Batf-directed CUT&Tag analysis revealed Batf binding signals in the open chromatin regions of the promoters for cytokine genes LIF and CLCF1, suggesting that Batf may potentially regulate these cytokine genes. This finding is consistent with previously reported inflammatory responses mediated by Batf2/Irf in macrophages (52). Based on the results from ATAC-seq, RNA-seq, and CUT&Tag analyses, it is suggested that AP-1 may potentially regulate the expression of numerous cytokines and participate in the regulation of the inflammatory response. This is consistent with previously reported findings that AP-1 plays a role in regulating the expression of inflammatory cytokines (47, 53–55).

Our study provides a comprehensive epigenetic landscape of chromatin accessibility in large yellow croaker in response to iridovirus infection. The significant changes observed in chromatin openness at immune-related genes and the identification of key transcription factors involved in the host immune response offer valuable insights into the host’s antiviral response mechanisms. These findings provide novel perspectives on the complex mechanisms underlying host antiviral responses and could contribute to the development of new antiviral strategies.

References

• 1

  WangDWuF-XSongD-DGaoH-QLiuX-ZCuiL-F. The China fishery statistical yearbook. 1st ed. Beijing: China Agriculture Press (2022).

  • Google Scholar

• 2

  WangX-WAoJ-QLiQ-GChenX-H. Quantitative detection of a marine fish iridovirus isolated from large yellow croaker, Pseudosciaena crocea, using a molecular beacon. J Virol Methods. (2006) 133:76–81. doi: 10.1016/j.jviromet.2005.10.025

  • CrossRef
  • Google Scholar

• 3

  ChenXHLinKBWangXW. Outbreaks of an iridovirus disease in maricultured large yellow croaker, Larimichthys crocea (Richardson), in China. J Fish Dis. (2003) 26:615–9. doi: 10.1046/j.1365-2761.2003.00494.x

  • CrossRef
  • Google Scholar

• 4

  QinPMunang’anduHMXuCXieJ. Megalocytivirus and other members of the family iridoviridae in finfish: A review of the etiology, epidemiology, diagnosis, prevention and control. Viruses. (2023) 15:1359. doi: 10.3390/v15061359

  • CrossRef
  • Google Scholar

• 5

  KuritaJNakajimaK. Megalocytiviruses. Viruses. (2012) 4:521–38. doi: 10.3390/v4040521

  • CrossRef
  • Google Scholar

• 6

  KwonWJChoiJCHongSKimYCJeongMGMinJGet al. Development of a high-dose vaccine formulation for prevention of megalocytivirus infection in rock bream (Oplegnathus fasciatus). Vaccine. (2020) 38:8107–15. doi: 10.1016/j.vaccine.2020.11.001

  • CrossRef
  • Google Scholar

• 7

  ZhangQGuiJ-F. Virus genomes and virus-host interactions in aquaculture animals. Sci China Life Sci. (2015) 58:156–69. doi: 10.1007/s11427-015-4802-y

  • CrossRef
  • Google Scholar

• 8

  GuoSWenZZhangXLiFCuiHLinXet al. Characterization of five caspase genes and their transcriptional changes in response to exogenous iridescent virus challenge in the whiteleg shrimp (Litopenaeus vannamei). Aquaculture. (2021) 534:736192. doi: 10.1016/j.aquaculture.2020.736192

  • CrossRef
  • Google Scholar

• 9

  BuenrostroJDGiresiPGZabaLCChangHYGreenleafWJ. Transposition of native chromatin for fast and sensitive epigenomic profiling of open chromatin, DNA-binding proteins and nucleosome position. Nat Methods. (2013) 10:1213–8. doi: 10.1038/nmeth.2688

  • CrossRef
  • Google Scholar

• 10

  KlemmSLShiponyZGreenleafWJ. Chromatin accessibility and the regulatory epigenome. Nat Rev Genet. (2019) 20:207–20. doi: 10.1038/s41576-018-0089-8

  • CrossRef
  • Google Scholar

• 11

  KongXWeiGChenNZhaoSShenYZhangJet al. Dynamic chromatin accessibility profiling reveals changes in host genome organization in response to baculovirus infection. PloS Pathog. (2020) 16:e1008633. doi: 10.1371/journal.ppat.1008633

  • CrossRef
  • Google Scholar

• 12

  HuSYangSLuYDengYLiLZhuJet al. Dynamics of the transcriptome and accessible chromatin landscapes during early goose ovarian development. Front Cell Dev Biol. (2020) 8:196. doi: 10.3389/fcell.2020.00196

  • CrossRef
  • Google Scholar

• 13

  SunYMiaoNSunT. Detect accessible chromatin using ATAC-sequencing, from principle to applications. Hereditas. (2019) 156:29. doi: 10.1186/s41065-019-0105-9

  • CrossRef
  • Google Scholar

• 14

  PritykinYvan der VeekenJPineARZhongYSahinMMazutisLet al. A unified atlas of CD8 T cell dysfunctional states in cancer and infection. Mol Cell. (2021) 81:2477–2493.e10. doi: 10.1016/j.molcel.2021.03.045

  • CrossRef
  • Google Scholar

• 15

  LiMLiJZhangYZhaiYChenYLinLet al. Integrated ATAC-seq and RNA-seq data analysis identifies transcription factors related to rice stripe virus infection in Oryza sativa. Mol Plant Pathol. (2024) 25:e13446. doi: 10.1111/mpp.13446

  • CrossRef
  • Google Scholar

• 16

  WuCZhaoXBabuVSYuanGWangWSuJet al. Distribution of mannose receptor in blunt snout bream (Megalobrama amblycephala) during the embryonic development and its immune response to the challenge of Aeromonas hydrophila. Fish Shellfish Immunol. (2018) 78:52–9. doi: 10.1016/j.fsi.2018.03.049

  • CrossRef
  • Google Scholar

• 17

  WuCDaiYYuanGSuJLiuX. Immunomodulatory effects and induction of apoptosis by different molecular weight chitosan oligosaccharides in head kidney macrophages from blunt snout bream (Megalobrama amblycephala). Front Immunol. (2019) 10:869. doi: 10.3389/fimmu.2019.00869

  • CrossRef
  • Google Scholar

• 18

  MaZ-YLiangJ-XLiW-SSunYWuC-SHuY-Zet al. Complement C3a enhances the phagocytic activity of B cells through C3aR in a fish. Front Immunol. (2022) 13:873982. doi: 10.3389/fimmu.2022.873982

  • CrossRef
  • Google Scholar

• 19

  HuR-DZhangWLiLZuoZ-QMaMMaJ-Fet al. Chromatin accessibility analysis identifies the transcription factor ETV5 as a suppressor of adipose tissue macrophage activation in obesity. Cell Death Dis. (2021) 12:1–11. doi: 10.1038/s41419-021-04308-0

  • CrossRef
  • Google Scholar

• 20

  CaoLLuoYGuoXLiuSLiSLiJet al. SAFA facilitates chromatin opening of immune genes through interacting with anti-viral host RNAs. PloS Pathog. (2022) 18:e1010599. doi: 10.1371/journal.ppat.1010599

  • CrossRef
  • Google Scholar

• 21

  CaoLLiuSLiYYangGLuoYLiSet al. The nuclear matrix protein SAFA surveils viral RNA and facilitates immunity by activating antiviral enhancers and super-enhancers. Cell Host Microbe. (2019) 26:369–384.e8. doi: 10.1016/j.chom.2019.08.010

  • CrossRef
  • Google Scholar

• 22

  WangXHWangKRNiePChenXHAoJQ. Establishment and characterization of a head kidney cell line from large yellow croaker Pseudosciaena crocea. J Fish Biol. (2014) 84:1551–61. doi: 10.1111/jfb.12386

  • CrossRef
  • Google Scholar

• 23

  CorcesMRTrevinoAEHamiltonEGGreensidePGSinnott-ArmstrongNAVesunaSet al. An improved ATAC-seq protocol reduces background and enables interrogation of frozen tissues. Nat Methods. (2017) 14:959–62. doi: 10.1038/nmeth.4396

  • CrossRef
  • Google Scholar

• 24

  ChuanJZhouAHaleLRHeMLiX. Atria: an ultra-fast and accurate trimmer for adapter and quality trimming. GigaByte. (2021) 2021:gigabyte31. doi: 10.46471/gigabyte.31

  • CrossRef
  • Google Scholar

• 25

  BrownJPirrungMMcCueLA. FQC Dashboard: integrates FastQC results into a web-based, interactive, and extensible FASTQ quality control tool. Bioinformatics. (2017) 33:3137–9. doi: 10.1093/bioinformatics/btx373

  • CrossRef
  • Google Scholar

• 26

  LangmeadBWilksCAntonescuVCharlesR. Scaling read aligners to hundreds of threads on general-purpose processors. Bioinformatics. (2019) 35:421–32. doi: 10.1093/bioinformatics/bty648

  • CrossRef
  • Google Scholar

• 27

  TarasovAVilellaAJCuppenENijmanIJPrinsP. Sambamba: fast processing of NGS alignment formats. Bioinformatics. (2015) 31:2032–4. doi: 10.1093/bioinformatics/btv098

  • CrossRef
  • Google Scholar

• 28

  RamírezFRyanDPGrüningBBhardwajVKilpertFRichterASet al. deepTools2: a next generation web server for deep-sequencing data analysis. Nucleic Acids Res. (2016) 44:W160–5. doi: 10.1093/nar/gkw257

  • CrossRef
  • Google Scholar

• 29

  LiuT. Use model-based analysis of chIP-seq (MACS) to analyze short reads generated by sequencing protein–DNA interactions in embryonic stem cells. In: KidderBL, editor. Stem cell transcriptional networks. Methods in molecular biology. Springer New York, New York, NY (2014). p. 81–95. doi: 10.1007/978-1-4939-0512-6_4

  • CrossRef
  • Google Scholar

• 30

  QuinlanARHallIM. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics. (2010) 26:841–2. doi: 10.1093/bioinformatics/btq033

  • CrossRef
  • Google Scholar

• 31

  YuGWangL-GHeQ-Y. ChIPseeker: an R/Bioconductor package for ChIP peak annotation, comparison and visualization. Bioinformatics. (2015) 31:2382–3. doi: 10.1093/bioinformatics/btv145

  • CrossRef
  • Google Scholar

• 32

  HeinzSBennerCSpannNBertolinoELinYCLasloPet al. Simple combinations of lineage-determining transcription factors prime cis-regulatory elements required for macrophage and B cell identities. Mol Cell. (2010) 38:576–89. doi: 10.1016/j.molcel.2010.05.004

  • CrossRef
  • Google Scholar

• 33

  Ross-InnesCSStarkRTeschendorffAEHolmesKAAliHRDunningMJet al. Differential oestrogen receptor binding is associated with clinical outcome in breast cancer. Nature. (2012) 481:389–93. doi: 10.1038/nature10730

  • CrossRef
  • Google Scholar

• 34

  RobinsonMDMcCarthyDJSmythGK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. (2010) 26:139–40. doi: 10.1093/bioinformatics/btp616

  • CrossRef
  • Google Scholar

• 35

  DobinADavisCASchlesingerFDrenkowJZaleskiCJhaSet al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. (2013) 29:15–21. doi: 10.1093/bioinformatics/bts635

  • CrossRef
  • Google Scholar

• 36

  PerteaMPerteaGMAntonescuCMChangT-CMendellJTSalzbergSL. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat Biotechnol. (2015) 33:290–5. doi: 10.1038/nbt.3122

  • CrossRef
  • Google Scholar

• 37

  LiaoYSmythGKShiW. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics. (2014) 30:923–30. doi: 10.1093/bioinformatics/btt656

  • CrossRef
  • Google Scholar

• 38

  LoveMIHuberWAndersS. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. (2014) 15:550. doi: 10.1186/s13059-014-0550-8

  • CrossRef
  • Google Scholar

• 39

  YuGWangL-GHanYHeQ-Y. clusterProfiler: an R package for comparing biological themes among gene clusters. Omics J Integr Biol. (2012) 16:284–7. doi: 10.1089/omi.2011.0118

  • CrossRef
  • Google Scholar

• 40

  KancherlaJYangYChaeHCorrada BravoH. Epiviz File Server: Query, transform and interactively explore data from indexed genomic files. Bioinformatics. (2020) 36:4682–90. doi: 10.1093/bioinformatics/btaa591

  • CrossRef
  • Google Scholar

• 41

  YuGWangL-GYanG-RHeQ-Y. DOSE: an R/Bioconductor package for disease ontology semantic and enrichment analysis. Bioinformatics. (2015) 31:608–9. doi: 10.1093/bioinformatics/btu684

  • CrossRef
  • Google Scholar

• 42

  MartinM. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet.journal. (2011) 17:10. doi: 10.14806/ej.17.1.200

  • CrossRef
  • Google Scholar

• 43

  ZhangQChengTJinSGuoYWuYLiuDet al. Genome-wide open chromatin regions and their effects on the regulation of silk protein genes in Bombyx mori. Sci Rep. (2017) 7:12919. doi: 10.1038/s41598-017-13186-6

  • CrossRef
  • Google Scholar

• 44

  KimAYoonDLimYRohHJKimSParkC-Iet al. Co-expression network analysis of spleen transcriptome in rock bream (Oplegnathus fasciatus) naturally infected with rock bream iridovirus (RBIV). Int J Mol Sci. (2020) 21:1707. doi: 10.3390/ijms21051707

  • CrossRef
  • Google Scholar

• 45

  BentsenMGoymannPSchultheisHKleeKPetrovaAWiegandtRet al. ATAC-seq footprinting unravels kinetics of transcription factor binding during zygotic genome activation. Nat Commun. (2020) 11:4267. doi: 10.1038/s41467-020-18035-1

  • CrossRef
  • Google Scholar

• 46

  SongDLianYZhangL. The potential of activator protein 1 (AP-1) in cancer targeted therapy. Front Immunol. (2023) 14:1224892. doi: 10.3389/fimmu.2023.1224892

  • CrossRef
  • Google Scholar

• 47

  QiaoYHeHJonssonPSinhaIZhaoCDahlman-WrightK. AP-1 is a key regulator of proinflammatory cytokine TNFα-mediated triple-negative breast cancer progression. J Biol Chem. (2016) 291:5068–79. doi: 10.1074/jbc.M115.702571

  • CrossRef
  • Google Scholar

• 48

  BiddieSCJohnSSaboPJThurmanREJohnsonTASchiltzRLet al. Transcription factor AP1 potentiates chromatin accessibility and glucocorticoid receptor binding. Mol Cell. (2011) 43:145–55. doi: 10.1016/j.molcel.2011.06.016

  • CrossRef
  • Google Scholar

• 49

  ZhaoX-FPanH-FYuanHZhangW-HLiX-PWangG-Het al. Increased serum interleukin 17 in patients with systemic lupus erythematosus. Mol Biol Rep. (2010) 37:81–5. doi: 10.1007/s11033-009-9533-3

  • CrossRef
  • Google Scholar

• 50

  KyriakisJM. Activation of the AP-1 transcription factor by inflammatory cytokines of the TNF family. Gene Expr. (1999) 7:217. Available at: https://pmc.ncbi.nlm.nih.gov/articles/PMC6174675/.

  • Google Scholar

• 51

  MonjePHernández-LosaJLyonsRJCastelloneMDGutkindJS. Regulation of the transcriptional activity of c-fos by ERK: A NOVEL ROLE FOR THE PROLYL ISOMERASE PIN1*. J Biol Chem. (2005) 280:35081–4. doi: 10.1074/jbc.C500353200

  • CrossRef
  • Google Scholar

• 52

  RoySGulerRPariharSPSchmeierSKaczkowskiBNishimuraHet al. Batf2/Irf1 induces inflammatory responses in classically activated macrophages, lipopolysaccharides, and mycobacterial infection. J Immunol Baltim Md 1950. (2015) 194:6035–44. doi: 10.4049/jimmunol.1402521

  • CrossRef
  • Google Scholar

• 53

  GuoZZhuYXiaoHDaiRYangWXueWet al. Integration of ATAC-seq and RNA-seq identifies MX1-mediated AP-1 transcriptional regulation as a therapeutic target for Down syndrome. Biol Res. (2023) 56:67. doi: 10.1186/s40659-023-00474-x

  • CrossRef
  • Google Scholar

• 54

  AtsavesVLeventakiVRassidakisGZClaretFX. AP-1 transcription factors as regulators of immune responses in cancer. Cancers. (2019) 11:1037. doi: 10.3390/cancers11071037

  • CrossRef
  • Google Scholar

• 55

  YuXWangYSongYGaoXDengH. AP-1 is a regulatory transcription factor of inflammaging in the murine kidney and liver. Aging Cell. (2023) 22:e13858. doi: 10.1111/acel.13858

  • CrossRef
  • Google Scholar