Characterization of aminopeptidase encoding gene anp-1 and its association with development in Caenorhabditis elegans

cDNA synthesis and multiple alignments

C. elegans hermaphrodites (N2 strain) and anp-1 deletion strain (RB804) were obtained from the Caenorhabditis Genetics Centre (CGC), University of Minnesota, USA. All C. elegans strains were cultivated on NGM plates at 20 °C (Stiernagle, 2006). Wild type worms were harvested with mixed stages and washed three times with cold 1×M9 buffer. After grinding in liquid nitrogen, total RNA was extracted from the clean worms using TRIZOL (Invitrogen, Carlsbad, CA, USA) as per manufacturer′s instructions. Then, cDNA was generated using the SuperScript™ III Reverse Transcriptase kit (Invitrogen, Carlsbad, CA, USA). There are two transcript isoforms of anp-1 (WBGene00011587), t07f10.1a and t07f10.1b, only the longest isoform (t07f10.1a) was considered in the present study. The specific forward and reverse primers designed for anp-1 were 5′-ATGGCCTCCGCCTACACTTATG-3′ and 5′-TTCAAAGGCTCTTCGAGATCTC-3′. A total of 30 cycles of polymerase chain reaction (PCR) were performed with the following protocol: 95 °C for 5 min, 94 °C for 30 s, 56 °C for 30 s, 72 °C for 3 min and 72 °C for 10 min. The PCR product was purified using the QIAquick® PCR Purification Kit (QIAGEN, Hilden, Germany), cloned into a pEASY-T1 plasmid (TransGen, Beijing, China) and named “pEASY-anp-1”. The recombinant vector was transformed into DH5α competent cells, from which clones were identified by colony PCR and sequenced for verification by company (The Beijing Genomics Institute (BGI), Shenzhen, China).

In this study, aminopeptidase sequence alignments were generated by CLUSTAL W (Larkin et al., 2007) and ESPript 3. A phylogenetic tree of the M1 aminopeptidase family (cd09601) was rendered using MEGA6 software and FigTree v1.4.3 (Dereeper et al., 2008). The transmembrane region and glycosylation residues were predicted using TMHMM Server v. 2.0 and NetOGlyc 4.0 Server, respectively. The amino acid sequence of ANP-1 was used for homology modeling to predict the three-dimensional (3D) model crystal structure. The 3D structure was checked by PROCHECK.

Aminopeptidase activity assay and production of polyclonal antibody against recombinant ANP-1

The conserved M1 aminopeptidase sequence of ANP-1 was cloned into the pET28a vector and expressed as a recombinant protein (rANP-1) in the Escherichia coli (E. coli) strain BL21 (DE3) cells induced by 1 mM isopropyl β-D-1-thiogalactopyranoside (IPTG) for 4 h at 37 °C. The forward and reverse primers designed for ANP-1 carrying restriction enzyme sites for BamHI or XhoI (in italics) were 5′-CGCGGATCCATCAGTTATCAATTGACAGTCA-3′ and 5′-CCGCTCGAGAGGATATCCCAGCTGTTTAG-3′, respectively. Briefly, bacterial lysate pellet was washed in Buffer A (2 M urea) and Buffer B (4 M urea) before being resuspended in Buffer C (8 M urea, 25 mM Na₂HPO₄, 50 mM Tris-HCl, 100 mM NaCl, 10% (v/v) glycerol, pH 7.5). Then the protein were purified using an Ni-NTA affinity column (GE Healthcare, Chicago, IL, USA). The affinity column was washed with 10 volumes of Buffer D (40 mM imidazole) and Buffer E (60 mM imidazole), and the fusion protein was eluted with Buffer F (8 M urea, 250 mM imidazole, at pH 7.5). The eluted protein was dialyzed against 1x PBS with a urea gradient (4 M, 2 M, 1 M and 0 M urea) and protein concentration was estimated by a BCA protein assay kit (Pierce, Appleton, WI, USA). The molecular weight of rANP-1 protein was confirmed by Coomassie-stained SDS-PAGE gels (12%) and western blotting. Then, the polyclonal antibody against rANP-1 was prepared according to the previous study (Munkhjargal et al., 2016).

Aminopeptidase activity and inhibitor sensitivity of rANP-1 were tested via standard assays using the leucine aminopeptidase (LAP) substrate L-leucine-p-nitroanilide (L-Leu-pNA; Sigma) as previously described (Roberts et al., 2013). Absorbance was measured at 405 nm using a Tecan Sunrise Microplate Reader.

Proteins extraction and identification

Mix-stage N2 worms were maintained and harvested using liquid culture method (Zhou et al., 2014). Worms were stored at 80 °C until use. The membrane proteins were extracted following the manufacturer’s instructions of Mem-PER™ Plus Membrane Protein Extraction Kit (Thermo Fisher, Waltham, MA, USA). Then, the crude glycoproteins of C. elegans were extracted using ConA according to the previous study (Redmond, Geldhof & Knox, 2004). The membrane proteins and glycoproteins were confirmed by SDS-PAGE and western blot assays.

Lifespan assay

The animals were synchronized by alkaline hypochlorite treatment (20 mL total volume with 1.25 mL of 5 M NaOH, 5 mL of NaClO and 13.75 mL of water). Synchronized eggs were incubated in sterilized M9 buffer for about 16–20 h at 20 °C (200 revolutions/min) until the L1 larval stage, and the L1 larvae were transferred to freshly prepared NGM plates for lifespan experiments.

At least 50 L1-stage worms were picked and transferred to 5-fluoro-20-deoxyuridine (FUDR)-NGM medium dishes and seeded with fresh bacteria for food. This was recorded as Day zero (Lezzerini et al., 2015). Surviving worms were removed every 24 h until all of the worms had died. A worm was considered dead when it did not respond to a gentle touch with a pick needle. Nematode lifespan was measured as the time till death for 50% of the worms. Experiments were repeated three times independently.

Scoring body size and determining brood size

Two percent agarose glass pads were made according to the method described in Liang, Xiong & Savage-Dunn (2013). Briefly, 50 μL NaN₃ (25 mM) was added to each agarose pad labeled with a unique color, about 50 live L4/ young adult worms or adult worms were anesthetized in this NaN₃ solution and then coverslips were placed on the top. Images were captured and measured using a Zeiss stereoscope (Discovery v.12) with supporting software.

One synchronized L1-stage (see above) nematode was randomly picked from a correspondingly labeled NGM plate and cultivated into a young adult larva (about 48 h). The young adults were then shifted daily to new dishes before laying their first egg (Yu et al., 2015). Adult nematodes were moved every 6 h until they stopped laying eggs. The total number of eggs laid by individual hermaphrodites was calculated. Ten worms per strain were used.

RNA extraction, library preparation, and illumina hiseq xten sequencing

Total RNA was extracted from the L4/young adult worms according to the approach described previously and genomic DNA was removed using DNase I (TAKARA). RNA-seq transcriptome librariy was prepared following TruSeq™ RNA sample preparation Kit from Illumina (San Diego, CA, USA) using 5 μg of total RNA. Libraries were size selected for cDNA target fragments of 200–300 bp on 2% Low Range Ultra Agarose followed by PCR amplified using Phusion DNA polymerase (NEB) for 15 cycles. After quantified by TBS380, paired-end RNA-seq sequencing library was sequenced with the Illumina HiSeq xten (2 × 150 bp read length) at Majorbio Bio-Pharm Technology Co., Ltd., Shanghai, China.

Read mapping, differential expression analysis and functional enrichment

The raw paired end reads were trimmed and quality controlled by SeqPrep (https://github.com/jstjohn/SeqPrep) and Sickle (https://github.com/najoshi/sickle) with default parameters. Then clean reads were separately aligned to reference genome with orientation mode using TopHat (http://ccb.jhu.edu/software/tophat/index.shtml, version2.0.0) software (Trapnell, Pachter & Salzberg, 2009). To identify DEGs between two different samples, the expression level of each transcript was calculated according to the fragments per kilobase of exon per million mapped reads (FRKM) method. RSEM (http://deweylab.biostat.wisc.edu/rsem/) was used to quantify gene abundances (Li & Dewey, 2011). R statistical package software EdgeR (Empirical analysis of Digital Gene Expression in R, http://www.bioconductor.org/packages/2.12/bioc/html/edgeR.html) was utilized for differential expression analysis (Robinson, McCarthy & Smyth, 2010). In addition, functional-enrichment analysis including GO and KEGG were performed to identify which DEGs were significantly enriched in GO terms and metabolic pathways at Bonferroni-corrected p-value ≤ 0.05 compared with the whole-transcriptome background. GO functional enrichment and KEGG pathway analysis were carried out by Goatools (https://github.com/tanghaibao/Goatools) and KOBAS (http://kobas.cbi.pku.edu.cn/home.do) (Xie et al., 2011).

RNA extraction and qRT-PCR analysis

A total of 500 ng RNA (from synchronized L4/ young adult worms) was used for qRT-PCR analysis and performed in triplicate with the TransScript II Green One-Step qRT-PCR SuperMix kit (TRANSGEN BIOTECH, Beijing, China) according to the manufacturer’s instructions. Thermal cycling and fluorescence detection were performed using a Roche LightCycler 96 System. The fold change of gene expression was calculated using the 2^−ΔΔCt method. The relative fold change of mRNA expression was normalized to the level of GAPDH expression. The primers designed for qRT-PCR are listed in Table S1.

Statistical analyses

All data were represented as mean ± SEM. Student’s t-test was used to analyze brood size, body size, and relative mRNA expression levels between N2 strain and RB804 strain. Median survival days and p-values were calculated using the log-rank method (Mantel-Cox). All statistical analyses were conducted using GraphPad Prism 5.0. In all cases, the two-tailed p-value was calculated and p ≤ 0.05 was considered significant.

Multiple alignments and 3D structure of ANP-1

To investigate the similarity of ANP-1 with other aminopeptidases, we cloned and sequenced the gene encoding C. elegans aminopeptidase ANP-1 (Fig. S1). The molecular weight (MW) of ANP-1 was about 110 kDa based on the translated property of coding sequence. Multiple alignment results showed that ANP-1 was a member of the M1 aminopeptidase family (Fig. S2). The signature functional motif “HEXXH-(X)₁₈-E” and the exopeptidase motif “GXMEN” were conserved across all members of the M1 aminopeptidase family. PROTEIN BLAST result manifested that ANP-1 possessed M1 aminopeptidase N-2 (M1-APN-2) domain as well as an adipocyte-derived leucine aminopeptidase (ERAP1) functional domain. Sequence analysis revealed that ANP-1 had 42%, 33% and 34% identity with H11, R03G8.4 and R03G8.6, respectively. M1 aminopeptidase proteins from different species were analyzed based on the neighbor-joining (NJ) method. The phylogenetic relationship analysis of amino acid sequence data for H11, ANP-1, and their homologs (R03G8.4 and R03G8.6) from C. elegans showed strong support for a single clade (Fig. 1).

Homology modeling of ANP-1 was conducted with SWISS-MODEL (Fig. 2A). The aminopeptidase N (APN) 5lg6.1.A was selected based on its identity (about 34%) to ANP-1. The 3D structure of ANP-1 consisted of 28 helices, 31 β-pleated sheets, and 12 random coil (Fig. S3). One zinc ion was bound in the active site of the protein and the zinc-binding residues (HIS₃₆₆, HIS₄₀₀, and GLU₄₁₉) are represented with sticks (Fig. 2B). The 3D structure of ANP-1 protein was checked with PROCHECK for its stability and was found to be stable.

The asparagine-modified residues (GalNAc O-glycosylation sites) prediction result implied that ANP-1 was a transmembrane-glycopeptide protein. The feasible transmembrane helix of ANP-1 was from 29aa (amino acid, aa) to 51aa based on Hidden Markov Model (Fig. S2).

Aminopeptidase activity of recombinant ANP-1

The putative aminopeptidase fragment was amplified and cloned into a pET-28a plasmid (Fig. 3A). The rANP-1 was expressed in its insoluble form in the E. coli strain BL21 (DE3) induced by IPTG. The rANP-1 protein was purified and refolded before assessing by SDS-PAGE (12%) and western blotting. The results implied that we obtained rANP-1 protein and, the molecular weight of rANP-1 was about 52.7 kDa, which was consistent with the predicted result (Fig. 3B).

Aminopeptidase activity of rANP-1, measured by standard aminopeptidase assay, showed that rANP-1 retained the enzyme catalytic function. The enzymatic activity of rANP-1 was increased at higher concentrations (Fig. 3C). The aminopeptidase activity of rANP-1 was inhibited to a similar degree by the inhibitors amastatin and bestatin and the reduction rates of amastatin and bestatin were about 55% and 46%, respectively (Fig. 3D). Both the activity and inhibition assays suggested that the rANP-1 was refolded appropriately into an active protein.

Proteins composition of membrane and ConA-binding fraction of C. elegans

The membrane and ConA-binding peptides were separated under reducing SDS-PAGE conditions, ranging in size from approximately 15 kDa – >170 kDa (Fig. 4A, lane 1) and 25 kDa – >170 kDa (Fig. 4B, lane 1). Several of membrane proteins reacted with anti-rANP-1 serum, and the MW of those proteins were between 70 kDa and 130 kDa (Fig. 4A, lane 2). Western blot strip of C. elegans ConA-binding proteins were incubated with anti-rANP-1 serum. However, the antibody recognized only one band (Fig. 4B, lane 2), the band has the size of 110 kDa, which was consistent with the predicted MW of ANP-1.

Deletion of anp-1 decreased the lifespan and affected the body size and brood size in C. elegans

As shown in picture 5, the median survival of N2 group and RB804 group were about 18 days and 16 days, respectively. Compared to the N2 strain, the anp-1 mutation significantly decreased the lifespan of RB804 strain, p < 0.001 (Fig. 5A). Meanwhile, the body size of L4/young adult worms and adult worms was observed and measured as 701.2 ± 5.905 μm (N2 group, n = 50), 642.5 ± 10.66 μm (RB804 group, n = 48), 1158 ± 11.73 μm (N2 group, n = 50), and 1145 ± 14.25 μm (RB804 group, n = 49), which indicated that the body size of RB804 worms was significantly lower (p < 0.001) compared with the normal worms (Fig. 5C). However, there was no significant difference in adult body size of C. elegans when anp-1 deleted (Fig. 5D). Representative images of animal body length were shown in Fig. S4. Furthermore, we found that knockout of anp-1 gene significantly decreased total number of eggs of C. elegans (Fig. 5B), the brood size of N2 group and RB804 group were 206.7 ± 8.82 (n = 10) and 150 ± 7.95 (n = 10), respectively.

DEGs analysis and functional annotation

To provide some insight into the effect on gene expression patterns when anp-1 deleted, we performed transcriptome sequencing and cluster analysis of DEGs by distance calculation algorithm between N2 strain and RB804 strain. All sequencing data were submitted in NCBI and can be accessed in the Short Read Archive (SRA) under accession number PRJNA548230. Among the genes identified, 184 genes were differentially expressed with 164 genes downregulated and 20 genes upregulated between N2 group and RB804 group, using a false discovery rate ≤0.05 and a fold change ≥2 as significant cutoffs (Table S2). As shown in Fig. 6, the heatmap profile indicated DEGs among different groups. Of the DEGs, 184 genes were annotated into 47 sub-categories belonging to the following three main GO categories: biological process (BP), cellular component (CC) and molecular function (MF) (Fig. 7A). The BP sub-categories, CC sub-categories, and MF sub-categories were 22, 14 and 11, respectively. In BP sub-categories, the vast majority was related to cellular process, single-organism process, metabolic process, biological regulation, regulation of BP, developmental process, multicellular organismal process, reproduction, reproductive process, CC organization or biogenesis, localization and response to stimulus. In CC sub-categories, genes for cell, cell part, organelle, macromolecular complex, organelle part, and membrane were the top six. Among the MF sub-categories, the majority of the GO terms were grouped into binding, catalytic activity. The detailed information for the annotations was in Table S3.

GO enrichment analysis placed 184 DEGs into 1,062 functional categories: meiotic spindle midzone assembly was the largest cluster, followed by meiotic cell cycle, negative regulation of microtubule depolymerization, and negative regulation of microtubule polymerization or depolymerization, etc. (Fig. 7B). Based on KEGG pathway mapping, 85 DEGs were classified into the following 45 KEGG functional categories (e.g., FoxO signaling pathway, Apoptosis-multiple species, Nitrogen metabolism, Arginine biosynthesis, Dorso-ventral axis formation, mTOR signaling pathway). A summary of the findings are presented in Fig. 7C and Table S4. According to further KEGG analysis of the 164 down-regulated genes and 20 up-regulated genes between N2 strain and RB804 strain, the top sub functional categories were FoxO signaling pathway genes, mTOR signaling pathway genes, and ABC transporters, respectively. The detailed annotation for the results were summarized in Fig. S5 and Table S5.

qRT-PCR assay

To validate expression patterns of DEGs and to determine the potential roles of the genes in regulation pathway, we confirmed their differential expression patterns between wild type and anp-1 knockout animals by qRT-PCR. Expression level of DEGs were consistent with those obtained by transcriptome sequencing results, confirming the accuracy of the transcriptome sequencing results reported in this study. Those genes, including ins-7, daf-18, lip-1, rskn-1, meg-3, mex-1, lin-22, vit-2, gna-2, hop-1, and klp-5 were down-regulated, however, unc-52, ftn-1, col-7, col-10, haf-4, haf-9, F49E2.5, dpy-1, and glt-7 were up-regulated (Fig. 8; Fig. S6). Based on the KEGG classified results, these down-regulated genes (lip-1 and rskn-1) were classified as MAPK signaling pathway sub-cluster. The longevity regulating pathway components, daf-18 and ins-7 were down-regulated in RB804 group compared with N2 group. Daf-18 and cyb-1 belonging to FoxO signaling pathway sub-cluster were down-regulated in RB804 group (Table S6).