C57BL/6 mice (10 weeks, male) used in this study were provided by the Experimental Animal Center of Nantong University. All experimental animals were housed in a specific pathogen-free barrier system under a 12-hour light/dark cycle at 23 ± 2°C and 60% humidity. All animal experiments were approved by the Animal Ethics Committee of Nantong University and followed the International guidelines (NIH) for the Care and Use of Laboratory Animals. The tester was qualified to have the Jiangsu Provincial Laboratory Animal Practitioner Post Certificate with a certificate No. 220212305.

Mice were caged together. The emergence of a vaginal plug indicates the appropriate time for mating. Fetuses of embryonic day 20 were identified. Depending on their birth time, 1-, 4-, 7-, 14-, 21-day, 1-, 3-, 6-, and 12-month-old mice were identified. The sciatic nerve crush injury in C57BL/6 mice was induced by applying constant pressure to the nerve with a No. 5 jeweler’s forceps for 15 s. Injured sites were marked using carbon power, and the marker center was used to distinguish proximal and distal nerve stumps.

Total RNA samples were extracted from fresh or frozen tissues using the well-established TRIzol method. Briefly, TRIzol (Cat#15596018, Invitrogen) was used to immerse tissue samples, and a tissue crusher was used to completely dissolve these samples. Chloroform and isopropanol were sequentially added to aid RNA extraction and crystallization, respectively. The RNA pellet was rinsed twice with 75% ethanol. Concentration and OD260/OD280 were determined using a NanoDrop 2000 spectrophotometer (Thermo Scientific, United States), and agarose gel electrophoresis was used to detect RNA integrity to ensure RNA quality for subsequent sequencing or qPCR.

The transcriptome sequencing and analysis were conducted by OE Biotech Co., Ltd. (Shanghai, China). RNA samples from mouse sciatic nerves were processed according to the manufacturer’s instructions (TruSeq^® Stranded mRNA Sample Preparation Guide, Part # 15031047 Rev. D. September 2012) to build the library. Briefly, 4 μg of DNase-treated total RNA was incubated with magnetic beads with Oligo (dT) to capture the mRNAs with poly-A tails. The mRNAs were broken into short fragments by addition of fragmentation buffer. With the fragmented mRNA as a template and random hexamers primers, the first strand of cDNA was synthesized. The second-strand cDNA was synthesized in the prepared second-strand synthesis system. The double-strand cDNA was purified using a kit. The purified double-strand cDNA was then subjected to end repair, A-tailing, and adapter ligation for final PCR amplification. The constructed library was quality-checked with an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, United States), and then sequenced using the Illumina HiSeq X Ten platform (Illumina, San Diego, CA, United States) to generate double-ended data of 125 or 150 bp.

Hisat2 (Kim et al., 2015) was used to align the Raw Reads with the specified reference genome, to obtain the position information of the reads on the genome. The number of reads matched to protein-coding genes in each sample was counted using the htseq-count software (Anders et al., 2015), and converted into protein-coding gene expression levels and FPKM values according to the formula. The DESeq (Love et al., 2014) software was used to standardize the number of counts for each gene in each sample and the expression level of each gene was estimated using the BaseMean value. Fold of difference was calculated. Negative binomial distribution test was adopted to test the significance of the difference in the number of reads. Finally, the differential expression fold and the difference significance test results (P-value) were obtained.

DNA removal and cDNA synthesis were performed using HiScript II Q RT SuperMix for qPCR (+gDNA wiper) kit (Cat#R223-01, Vazyme, China). Reactions were as follows: first total RNA (1 μg), 4 × g DNA Wiper Mix (4 μL), and DEPC H2O were mixed into 16 μL, at 42°C for 2 min to remove DNA; 4 μL of 5 × HiScript II qRT SuperMix II was added at 50°C for 1 min and then at 85°C for 5 s to synthesize cDNA. Quantitative PCR reactions were performed on cDNA using ChamQ Universal SYBR qPCR Master Mix (Cat#Q711-02, Vazyme, China) in StepOnePlus™ Real-Time PCR System sequence detection system (Thermo Fisher Scientific). The 2^–ΔΔCt method was used to calculate relative changes of gene expression. Actb was used as the reference gene for PCR normalization. The primers are as follows:

Considering gene expression changes during development, injury, and regeneration, a P-value less than 0.05 indicates no differential expression and a P-value higher than 0.05 indicates differential expression. Three modules described in Figure 3 were defined based on the rule that “the change trend of target gene expression during development is consistent with that during regeneration but opposite to that during injury.” There are two possibilities for a gene to be judged as module 1: the first is that all three processes follow the above rule, and the second is that two processes meet the above rule with the third process unchanged. Genes with expression changes that do not conform to the above rule in ≥2 processes are judged as module 2. Genes with no expression changes in ≥2 processes were classified as module 3. Therefore, module 1 indicates basically consistent, module 2 indicates discrepancy, and module 3 indicates no change.

Fresh tissue samples were immersed in protein lysis buffer containing protease inhibitors and crushed using a crusher until completely dissolved. Lysate supernatant was mixed with loading buffer at a ratio of 4:1, and the mixture was boiled at 95°C for 10 min to obtain whole protein supernatant. Protein samples were quantified by spectrometer to ensure a total protein of 20 μg per lane. Protein samples separated by SDS-PAGE gel electrophoresis were transferred to wetted PVDF membranes (Millipore, Bedford, MA, United States) using the common Western blot method. To measure the amount of acetylated tissue, immunoblotting assay with anti-acetyllysine mouse mAb (Cat# PTM-101, PTM BIO, China) as primary antibody and Goat Anti-Mouse IgG H&L (Cat# ab205719, Abcam, United States) as secondary antibody was performed. The immunoblots were visualized using the enhanced chemiluminescence kit (Beyotime, Shanghai, China). Immunoblot intensities were quantified by ImageJ software.¹ The total gray value of all detectable bands in each lane was used for pan-acetylation quantification in Figure 4. Ponceau S staining was used to visualize total protein across the membrane.

Two databases were used to analyze the potential regulatory networks of Hdac1 and Kat2b. RRUST database (Han et al., 2018) was used to analyze the interaction of Hdac1 and Kat2b as transcription factors with the target, and the Gene Ontology biological process terms of the target were also provided by this database. The protein–protein interactions and interaction scores of Hdac1 and Kat2b with the target were downloaded from the STRING database (Szklarczyk et al., 2017). For interaction scores, Score > 200 was taken as cut-off. Venny was used to identify overlaps between downstream molecules (also included in Figure 5). DAVID (Huang et al., 2009) was used to analyze the Gene Ontology biological process terms of overlapping molecules. With target as node and score as edge, Cytoscape software² was used to generate interaction diagrams.

All statistical analyses were performed using SPSS 20.0 (IBM, United States) and experimental data were presented as SD ± SEM. SPSS Pearson analysis was used to calculate R and P-values. Intergroup comparison was done using one-way analysis of variance and P-values were thus calculated.

We prepared relevant animal models and conducted deep sequencing to systematically analyze the expression changes of acetyltransferases and deacetylases during sciatic nerve development and regeneration. For the developmental process, sampling was done at intensive time points, including embryonic day 20, postnatal days 1, 7, 14, and 21, postnatal months 1, 3, 6, 9, and 12, throughout the embryonic, postnatal, juvenile, and adult stages. For the regeneration process, we prepared the animal model of sciatic nerve clamping and set five observation time points, including postoperative days 1, 3, 7, 14, and 21, for proximal (proximal nerve stump, p) and distal (distal nerve stump, d) sampling according to the location of the injury. Total RNA was extracted to generate RNA-seq data. Each sample contained 40–50 million raw reads, and the Q30 value of all samples was >91.4%, suggesting that the sequencing depth and data are sufficiently available for subsequent analysis. The expression of all detected genes in developmental and injury samples is listed in Supplementary Files 1, 2.

In this study, we identified the coding genes of 20 deacetylases and 13 acetyltransferases and analyzed their expression changes using the Cluster heatmap. During development, deacetylases could be clearly divided into two categories, one containing most of the Hdac family members that were gradually down-regulated and the other containing most of the Sirt family members that were gradually up-regulated (Figure 1A). There was no uniform expression characteristics of acetyltransferases. The expression of Kat2b/5/8 was gradually up-regulated, the expression of Atat1, Esco2, Hat2, Kat2a, and Kat7 was gradually down-regulated, and the expression of other genes first increased and then decreased (Figure 1B). Overall, almost all genes expressed very similarly at E20 and 1d, and then their expression changed in a constant state until the greatest amount of change was achieved and stabilized at 6–9 months after birth (Figures 1A,B). To ensure the reliability of the sequencing results, we randomly selected five of the above genes to perform reverse transcription quantitative polymerase chain reaction (RT-qPCR) assay. Correlation coefficients were used for evaluating the consistency between experimental and sequencing results, and a value of R² = 0.86 indicated a high correlation between the two (Figure 1C).

Expression of acetyltransferases and deacetylases showed a complex variation with no discerning trend during sciatic nerve injury and regeneration. Unsurprisingly, the gene expression patterns of proximal and distal nerve stumps were not identical, which is consistent with the different physiological changes of proximal and distal nerve stumps (Min et al., 2021). However, there were clear and noticeable gene expression transitions for both acetylases and deacetylases at two periods, N-1d (Figures 2A–D; green boxes) and C7d-C14d (Figures 2A–D; blue boxes), corresponding to the injury period and regeneration periods. This implicates that acetylation shows more dramatic changes at these two periods than the other time periods. Interestingly, during the injury time period (N-1d) of the distal nerve stump samples, the expression pattern of Sirt genes is as follows: most of the Sirt family members were down-regulated, while most of the Hdac family were up-regulated (Figure 2A). This trend, which is opposite to that during development (Figures 1A vs. 2A), suggests a possible link between the expression of acetylated regulators during development, injury, and regeneration. We performed qPCR to verify this part of the sequencing results, and the two also showed high correlation (Figure 2E).

Based on the above clues, we hypothesized that there are common rules for the changes in the expression of acetylated regulators during sciatic nerve development, injury, and regeneration. We independently defined three gene expression modules for sorting acetyltransferase and deacetylase encoding genes during three time periods: 1d-6m (development), N-C1d (injury), and C7d-C14d (regeneration). According to the rule that “the target gene shows the same changes in the development and regeneration process but an opposite trend in the injury process,” module 1 indicates the expression of the target gene basically follows this rule, module 2 indicates the expression of the target gene does not follow this rule, and module 3 indicates no expression change (more details in section “Materials and Methods”). Modular classifications were conducted four times according to the type of each enzyme and the location of each sample. Genes and their expression changes contained in each module were displayed on a heat map (Figures 3A,D,G,J). The proportion of differentially and non-differentially expressed genes was used to assess whether acetylation and deacetylation were activated during this time period (Figures 3B,E,H,K). Deacetylation was found to be more active than acetylation during development (Figures 3B,H vs. E,K). Acetylation was similar between distal and proximal nerve stumps during injury and regeneration (Figures 3E,K), but deacetylation was more active distally than proximal (Figures 3B,H). Based on the number and proportion of genes in each module, we found that regardless of acetylation or deacetylation, module 1 had the largest proportion (Figures 3C,F,I,L). Our results reveal that the genes in module 1 are decisive, that is, most of acetylation regulators have the same changing trend during development and regeneration and the opposite trend during injury, which also implies that protein acetylation is involved in the above three processes.

In addition to expression analysis of acetylation regulators, we also tested changes in overall acetylation levels of the sciatic nerve during development, injury, and regeneration using pan-acetylation antibodies. Total acetylation levels in nerve tissue increased substantially from early development (1 day) to adulthood (6 months), which were approximately threefold higher at 6 months than at 1 day (Figure 4A). Subsequently, we constructed a mouse sciatic nerve clamp model, and Western blot assay was done at 1, 7, and 14 days to analyze the level of pan-acetylation. As expected, the total acetylation level during sciatic nerve injury (N vs. 1d) were ∼2.5-fold down-regulated as opposed to that during development (Figure 4B). During post-injury regeneration (7d vs. 14d), the total acetylation level in sciatic nerve tissue was approximately 1.6-fold up-regulated, which was consistent with that during development but opposite to that during injury (Figure 4C). Our results suggest that changes in overall acetylation during sciatic nerve development, injury, and regeneration are consistent with the regularity of acetylation regulators and may imply that the upregulation of acetylation during development is primarily driven by downregulated deacetylases and upregulated acetyltransferase.

We sought to find which regulators might act as major genes associated with development, injury, and regeneration. The Venny tool was used to intersect differentially expressed genes during development, injury, and regeneration. Ultimately, Hdac1/2, Kat2b, and Esco2 emerged as differentially expressed genes at any part of any process (Figures 5A–D). Upon examination, Hdac1/2 was down-regulated and Kat2b was up-regulated during development, consistent with the trend of increased pan-acetylation during development.

To uncover the mechanisms underlying key acetylation regulators associated with sciatic nerve development, injury, and regeneration, we analyzed their potential downstream molecules represented by Hdac1 and Kat2b. Acetyltransferases and deacetylases act by regulating both histone acetylation and non-histone acetylation (Marmorstein and Zhou, 2014; Sabari et al., 2017). Therefore, we analyzed the regulatory molecules of Hdac1 and Kat2b from the perspectives of transcriptional activation and protein interaction. A total of 28 and 5 TF-target interaction pairs for Hdac1 and Kat2b were identified respectively using the Transcriptional Regulatory Relationships Unraveled by Sentence-based Text Mining (TRRUST). As transcription factors, the downstream molecules of Hdac1 and Kat2b did not overlap. Functional analysis of Hdac1-targeted molecules indicated that they could be clustered for the regulation of cell proliferation and cell cycle. In contrast, Kat2b-targeted molecules could be aggregated in the cell differentiation pathway (Figure 6A). These pathways were all related to cell cycle. We also used STRING to obtain 3,548 protein–protein interaction pairs for Hdac1 and 2,967 for Kat2b. We found that these two proteins had a general interaction with other acetylation regulators, suggesting that a complex regulatory relationship, co-regulated by multiple enzymes, exists in the intracellular acetylation. Hdac1 showed stronger interaction ability and regulated more molecules than Kat2b (Figure 6B). Unlike transcription factor activity, most of the downstream molecules of Hdac1 and Kat2b were overlapped for protein interaction, suggesting that Hdac1 and Kat2b may work by regulating the same molecules (Figure 6C). GO analysis of these 2,291 molecules indicated that, unsurprisingly, some of the most prominent pathways were mostly clustered in the transcription regulation pathway regulated by histone acetylation. However, for non-histone regulation, most of these molecules can be aggregated in cell cycle-related pathways such as regulation of cell proliferation, cell cycle, cell apoptosis, and cell differentiation (Figure 6D). We focused on the interaction between the two proteins and proteins related to the cell proliferation regulation pathway. The interaction of Hdac1 with this pathway was slightly stronger than that of Kat2b. Among the proteins related to the cell proliferation regulation pathway, Dmnt1, Trp53, Smad4, Mbd2, E2f4, Ctnnb1 interacted more strongly with Hdac1, and Notch1, Trp53, Rbpj, Sirt1, Tal1, Nr3c1 interacted more strongly with Kat2b (Figure 6E). Importantly, these molecules have been clearly reported to have acetylation modifications that affect their own functions (Huang et al., 2014; Kong et al., 2020; Jia et al., 2021). Notch1, Cdk4, Sirt1, and Nr3c1 are confirmed to be closely related to peripheral nerve regeneration or SCs physiological functions (Woodhoo et al., 2009; Viader et al., 2011; Lerch et al., 2017; Romeo-Guitart et al., 2019).

Changes in the expression of downstream molecules shared by Hdac1 and Kat2b in the cell proliferation pathway were analyzed during development, injury, and regeneration. These target genes changed most dramatically during development, followed by injury and milder regeneration (Figures 7A–C). Regardless of injury or regeneration, differentially expressed targets are similar in proximal and distal nerve stumps, consistent with the expression signature of acetylation regulators (Figures 3B,E,H,K, 7B,C). Notably, IL6 is significantly differentially expressed during the injury, suggesting that acetylation may be one of the mechanisms of inflammatory activation after peripheral nerve injury (Figure 7B). The expression changes of most target genes during development, injury, and regeneration are consistent with “development and regeneration are consistent and opposite to injury,” among which representative targets are Foxo1, Cdkn1b, Nr3c1, Jup, Stat6, Jund, and Hck.

Our data suggest that key acetylation regulators such as Hdac1 and Kat2b may mediate morphological changes that are similar or opposite during sciatic nerve development, injury, and regeneration through regulating the acetylation levels of cell cycle pathway-related members (Figure 8).