Analysis of RNA-Seq data reveals involvement of JAK/STAT signalling during leg regeneration in the cricket Gryllus bimaculatus

Regeneration requires the recognition of tissue loss and the subsequent restoration of the relevant structure. Several research areas have clarified the molecular bases of tissue regeneration in model organisms (Agata and Inoue, 2012; Tanaka and Reddien, 2011). Key molecules for the early regeneration process include apoptotic factors and proteinases, such as matrix metalloproteinases (MMPs) that regulate cell motility or morphology via degradation of the cell surface and cell-cell contact (Chera et al., 2011; Maki et al., 2010; Repiso et al., 2011).

Pluripotent stem cells or dedifferentiated multipotent cells proliferate in response to injury to form a blastema at the tip of the amputated stump. A blastema is a population of pluripotent proliferating cells that can restore the lost parts of tissues. Signalling molecules, including TGFβ, Shh, Wnt, FGF, EGF and IGF, are expressed in the epidermis of the wound, and in some cases in blastema cells, to regulate the proliferation and differentiation of blastema cells via the MAPK cascade (Chera et al., 2011; Petersen and Reddien, 2009; Repiso et al., 2011; Satoh et al., 2011; Tasaki et al., 2011a; Tasaki et al., 2011b; Wang et al., 2009; Yazawa et al., 2009; Yoshinari and Kawakami, 2011). The activities of MAPKs are negatively regulated by MAPK phosphatases (Repiso et al., 2011; Tasaki et al., 2011a; Tasaki et al., 2011b). In parallel with EGF signalling, the Salvador/Hippo/Warts and JAK/STAT signalling pathways regulate the proliferation of stem cells (reviewed by Jiang and Edgar, 2011).

The two-spotted cricket Gryllus bimaculatus is an emerging model organism for regeneration research (Mito and Noji, 2009). When we amputate the leg of the cricket nymph, a blastema forms at the distal tip of the leg stump to restore the lost part of the leg (Mito et al., 2002). The early blastema displays cell proliferation and activation of several signalling pathways, such as Wg/Wnt and EGFR, at 2 days postamputation (dpa) (Mito et al., 2002; Nakamura et al., 2008b). However, the molecular mechanisms of blastema formation remain to be identified.

In the present study, we analysed the global gene expression profile during the early regeneration process to identify the molecules involved in blastema formation. We found that the expression levels of components of the JAK/STAT signalling pathway were upregulated in blastemas, and that JAK/STAT signalling was essential for leg regeneration. We conclude that comparative transcriptome analysis could be useful to clarify the molecular mechanisms of blastema formation.

G. bimaculatus used in this study were reared under standard conditions (Mito and Noji, 2009).

The regenerating tibiae or blastemas of control or RNAi-treated nymphs (n=9-15) were pooled into single tubes and total RNA was extracted using the RNAqueous-Micro Kit (Ambion). Each pooled RNA was divided into two samples and each half reverse transcribed to prepare cDNA. Each cDNA was divided into three samples and used for q-PCR. Relative transcript ratios in the q-PCR study were calculated from experiments performed in triplicate and are shown as mean ± s.d. q-PCR experiments were repeated twice for confirmation. The total number of nymphs used to extract RNA is indicated by n. The housekeeping gene Gryllus β-actin was used as an internal control (Bando et al., 2011a; Bando et al., 2009).

Distal parts of regenerating legs at 0 and 24 hours postamputation (hpa), designated as normal legs (NLs) and regenerating legs (RLs), respectively, were separately collected. The cDNA libraries constructed from poly(A)⁺ RNAs were sequenced using the Illumina Genome Analyzer IIx next-generation sequencer by a 50 bp single-ended process. Raw sequence data from NLs and RLs have been deposited in the DNA Data Bank of Japan (DDBJ) sequence read archive (accession numbers of NL and RL are DRR001985 and DRR001986, respectively).

To construct the assembled transcripts, all of the reads obtained from the NL and RL samples were assembled using the Velvet (EMBL) and Oases (EMBL) de novo transcriptome assembly software packages. Each read was mapped to the assembled transcripts using the Bowtie short-read aligner (Johns Hopkins University). Transcript expression levels were estimated in order to calculate the reads per kb of exons per million mapped reads (RPKM) values. Assembled transcripts that were upregulated more than twofold in RLs compared with NLs based on the RPKM values were annotated with the BLASTX program against the NCBI nonredundant protein sequence database, with an E-value cut-off of 0.001.

Functional annotation and Gene Ontology (GO) term assignment were performed with the GO-mapping algorithm of the Blast2GO program (Conesa et al., 2005) and the InterProScan sequence search tool (EMBL) with default settings. The GO annotations obtained by these two programs were merged. The GO-Slim view was used before generating multilevel GO graphs.

Gryllus homologues of domeless, hopscotch, Stat, Socs2, Socs5, Socs7 and su(var)2-10 [Gb’dome1/2 (Gb’dome), Gb’hop1/2 (Gb’hop), Gb’Stat1/2 (Gb’Stat), Gb’Socs2, Gb’Socs5, Gb’Socs7 and Gb’su(var)2-10, respectively] were cloned by PCR using cDNA from the whole body of cricket nymphs at third instar. The cDNA sequences have been deposited in GenBank (accession numbers AB715378-AB715379 and AB761612-AB761618). Preparation of dsRNA and RNAi were performed as previously described (Bando et al., 2009). All RNAi data are given as mean ± s.d. and represent three to five independent experiments. The total numbers of RNAi-treated nymphs of multiple independent experiments are indicated by n.

The early blastema is formed within 2 dpa during the leg regeneration process, as revealed previously by cell proliferation assays and the expression of Gb’Egfr (Mito et al., 2002; Nakamura et al., 2008b). We studied the temporal changes in the mRNA levels of Gb’Egfr and Gb’cyclinE by q-PCR (supplementary material Fig. S1C). The expression of Gb’Egfr and Gb’cyclinE doubled from 12 to 132 hpa and from 24 to 132 hpa, respectively, compared with the levels at 0 hpa (supplementary material Fig. S1C), suggesting that EGF-mediated blastema formation starts from 24 hpa.

To uncover the molecular basis of blastema formation, we performed a transcriptome analysis of the early blastema in G. bimaculatus without a reference genome (supplementary material Fig. S1D). Early blastemas of RLs at 24 hpa and the corresponding leg part of NLs, in which the distal portion of the leg was cut off just prior to experiments, were collected from sixth instar nymphs (supplementary material Fig. S1A,B). The cDNAs synthesized using RNA samples from RLs and NLs were sequenced in two independent lanes using a next-generation sequencer.

All sequence reads from both RLs and NLs were assembled into 55,043 transcripts using transcriptome assembly software packages (supplementary material Fig. S1, Table S1). The average length of the assembled transcripts was 261 bp (supplementary material Fig. S2A, Table S1). Each read obtained from the RLs and NLs was mapped to the assembled transcripts to calculate the RPKM value of each assembled transcript (see Materials and methods). A total of 466 transcripts were only expressed in RLs, and the expression of 11,492 transcripts was upregulated more than twofold in RLs (supplementary material Fig. S1D, Table S1) based on the comparison of RPKM values between RL and NL. Because the transcripts that were upregulated or only expressed in RLs could be involved in blastema formation, we focused on these 11,958 transcripts for further annotation.

We analysed the 11,958 transcripts using the Blast2GO program (supplementary material Fig. S1D). In total, 2724 transcripts had significant hits with known genes annotated by BLASTX (supplementary material Table S1). Gryllus transcripts showed high homologies to the orthologous genes of Pediculus humorous and Tribolium castaneum (supplementary material Fig. S2B). Only 23% of the transcripts were annotated because the length of the assembled transcript obtained in this work was short and the genome or EST data are currently insufficient in G. bimaculatus.

To obtain a general overview of the functions of the transcripts that were upregulated during regeneration, GO terms were assigned. Of the 2724 transcripts, 2344 were annotated with at least one GO term (supplementary material Table S1). Frequently encountered GO terms for biological process, molecular function and cellular components (supplementary material Fig. S3) perhaps relate to the fact that signalling processes transduce extracellular signals received on the cell membrane into the nucleus via protein-protein interactions and protein modifications.

We searched the assembled transcripts to identify Gryllus homologues of components known to be involved in regeneration in other animals (Agata and Inoue, 2012; Nakamura et al., 2008a). The assembled transcripts encoding candidate genes are listed in Table 1 and supplementary material Tables S2-S6. A high RL/NL RPKM value was given if the candidate gene was encoded by several assembled transcripts.

The expression of apoptotic genes, MMPs, piwi and AGO3 was upregulated in RLs (supplementary material Table S2), similar to regeneration in other organisms (Chera et al., 2011; Li et al., 2011; Maki et al., 2010; Reddien et al., 2005). Various signalling components, including Wg/Wnt, Hh/Shh, EGF, Dpp/TGF, VEGF, Insulin/IGF, Notch, Toll (supplementary material Table S3), Dachsous/Fat (supplementary material Table S4) and JAK/STAT (Table 1) were also upregulated in RLs (Agata and Inoue, 2012; Bando et al., 2011a; Bando et al., 2009; Bando et al., 2011b; Nakamura et al., 2008a; Nakamura et al., 2008b). Supplementary material Tables S5 and S6 list transcription and epigenetic factors, respectively, that were upregulated in RLs, which could be key factors for cell fate determination and reprogramming (Maki et al., 2010; McClure and Schubiger, 2008; Meissner, 2010; Onder et al., 2012; Rao et al., 2009; Repiso et al., 2011; Stewart et al., 2009).

Toll, Dpp/TGFβ and JAK/STAT were the top three signalling pathways upregulated in RLs (Table 1; supplementary material Table S3). The evolutionarily conserved JAK/STAT signalling pathway has roles in embryogenesis, immunity and cell proliferation (Arbouzova and Zeidler, 2006) and is involved in the proliferation and maintenance of stem cells in Drosophila (Jiang and Edgar, 2011; Karpowicz et al., 2010; Ren et al., 2010; Shaw et al., 2010; Staley and Irvine, 2010). Thus, activation of JAK/STAT signalling might be required for blastema formation.

To confirm the transcriptome data, we estimated the transcription levels of Gryllus homologues of the interleukin receptor Domeless (Gb’dome), the Janus kinase Hopscotch (Gb’hop) and the transcription factor STAT (Gb’Stat) during regeneration by q-PCR. Gb’dome, Gb’hop and Gb’Stat were upregulated at 12-24 hpa (Fig. 1B). The relative ratios of these genes at 24 hpa were 2.1±0.2, 1.9±0.2 and 2.1±0.2, respectively, suggesting that our transcriptome analysis was adequate.

Two independent regions of each of the three genes were cloned (supplementary material Fig. S4A). RNAi was performed (Fig. 1A) to determine whether JAK/STAT signalling is involved in early blastema formation. The relative ratios of endogenous Gb’dome, Gb’hop and Gb’Stat transcripts were lower (P<0.01, Student’s t-test) in the respective RNAi crickets compared with the control at 5 dpa as assessed by q-PCR, indicating that RNAi had occurred (supplementary material Fig. S5C). We compared the phenotypes obtained with two dsRNAs corresponding to different regions of each gene. We observed a similar defect during leg regeneration with RNAi against Gb’dome, Gb’hop or Gb’Stat (Fig. 1C; supplementary material Fig. S5A,B). This finding indicates that the phenotypes obtained by RNAi were not off-target effects.

In control experiments, when the metathoracic tibia of the third instar nymph was amputated, a blastema formed at the distal end of the amputated leg during the third to fourth instar stages. In the fifth instar, miniaturised forms of the tibiae and tarsi were restored. In the sixth instar and adult stages, the amputated legs restored the missing portion to regain a nearly normal appearance (Fig. 1A). Gb’dome^RNAi, Gb’hop^RNAi and Gb’Stat^RNAi nymphs were viable and showed defects during leg regeneration. In the fourth instar, Gb’dome^RNAi, Gb’hop^RNAi and Gb’Stat^RNAi nymphs showed no obvious defects. However, tissue regeneration of the lost parts of the legs did not occur in the fifth and sixth instar stages (Fig. 1A; supplementary material Fig. S5A).

RNAi phenotypes were categorised into four classes at sixth instar (Fig. 1C). The major phenotypes of each RNAi experiment are shown in Fig. 1A. In the Gb’dome^RNAi crickets, distal structures, including tibial spurs, were regenerated. However, the lost parts of the tarsi were not regenerated (Fig. 1A, asterisk). In Gb’hop^RNAi crickets, the distal portion of the tibia, shown by tibial spurs, was regenerated. The lost tarsus was regenerated as a nonsegmental structure in miniature (Fig. 1A, arrowheads) in the sixth instar and adult stages. In Gb’Stat^RNAi crickets, lost parts of the legs were not regenerated (Fig. 1A). Given that Dome transduces signals to other pathways, that STAT is phosphorylated by other kinases (Arbouzova and Zeidler, 2006) and that the RNAi efficiencies against Gb’dome, Gb’hop and Gb’Stat differed, it is conceivable that the phenotype of Gb’hop^RNAi was milder than those of Gb’dome^RNAi and Gb’Stat^RNAi.

We compared the lengths of the regenerated and normal tibiae at sixth instar. In control regenerations, when tibiae were amputated at distal (70%) positions, the lengths of the regenerated tibiae were 101±5% (n=22) of unamputated control tibia lengths (supplementary material Fig. S5A,B). In the Gb’dome^RNAi, Gb’hop^RNAi and Gb’Stat^RNAi nymphs, the lengths of the regenerating tibiae were shortened by distal amputation; the regenerated tibiae were 66±7% (n=23), 80±4% (n=22) and 64±8% (n=19) of the control tibia lengths, respectively (P<0.01; supplementary material Fig. S5A,B). The lengths of the regenerated tibiae were similar to that of the amputated tibia at third instar; thus, Gb’dome and Gb’Stat appear to regulate only blastema cell proliferation but not allometric growth of the leg stump.

To analyse the function of JAK/STAT signalling further, we focused on the function of Gb’Stat. The Gb’Stat^RNAi nymphs showed normal body morphologies and weights (111.4±15.6 mg, n=4) compared with control nymphs (109.3±16.3 mg, n=11, P>0.05; Fig. 2A,B), except for the morphologies of the regenerated legs (Fig. 2D). Because the Gb’Stat^RNAi cricket showed defects in leg regeneration but not in the morphology of other tissues, we sought to clarify whether Gb’Stat promotes blastema cell proliferation. We attempted to detect proliferating cells using EdU, but we could not remove cuticles without damaging the regenerating Gb’Stat^RNAi legs, as these were more fragile than control legs. EdU is incorporated into S-phase cells, when cyclin E is specifically expressed, so we determined the ratio of proliferating cells by Gb’cyclinE expression using q-PCR as an alternative (Bando et al., 2011a). The relative ratio of Gb’cyclinE mRNA in the blastemas at 5 dpa was decreased in Gb’Stat^RNAi nymphs (56±5%, n=6, P<0.01; Fig. 2C) compared with control nymphs (100±10%, n=7). This suggests that Gb’Stat could promote blastema cell proliferation.

We performed dual RNAi experiments against Gb’Stat with Gb’fat (Gb’ft) or Gb’warts (Gb’wts). We previously showed that blastemas are enlarged by hyperproliferation of blastema cells in Gb’ft^RNAi and Gb’wts^RNAi nymphs at fourth instar (Fig. 2E, asterisks) (Bando et al., 2009). Hyperproliferation in the blastema caused by Gb’ft^RNAi or Gb’wts^RNAi was suppressed by RNAi against Gb’Stat (Fig. 2E). The regeneration-defective phenotype caused by Gb’Stat^RNAi was dominant to the abnormal regeneration phenotypes caused by Gb’ft^RNAi or Gb’wts^RNAi at the sixth instar and adult stages (Fig. 2F,G). Therefore, Gb’Stat appears to promote blastema cell proliferation via expression of Gb’cyclinE during Gryllus leg regeneration.

The transcription of Gb’dome, Gb’hop and Gb’Stat was not uniformly regulated by Dachsous/Fat signalling. The transcription of Gb’Stat was decreased in Gb’ft^RNAi and Gb’wts^RNAi nymphs, whereas that of Gb’dome and Gb’hop was increased (supplementary material Fig. S6). This indicates that the JAK/STAT and Dachsous/Fat signalling pathways may not be epistatic.

Our comparative transcriptome analysis indicates that several negative regulators of JAK/STAT signalling were also upregulated in RLs (Table 1). We performed RNAi against Gb’Socs2, Gb’Socs5, Gb’Socs7 and Gb’su(var)2-10 (Fig. 3A; supplementary material Fig. S4A,B). In all experiments, the amputated surfaces of legs were covered by cuticle at fourth instar, and miniaturised forms of the tarsi and tarsal claws were restored at fifth instar. In the sixth instar and adult stages, the missing portions of the amputated legs were restored, and a nearly normal appearance was regained (Fig. 3A). However, the regenerating tibiae of Gb’Socs2^RNAi (110±8%) and Gb’Socs5^RNAi (108±5%) crickets were significantly lengthened, and those of Gb’su(var)2-10^RNAi (82±12%) crickets were significantly shortened (Fig. 3A; supplementary material Fig. S5A,B).

Because SOCS negatively regulates JAK/STAT signalling activity, we sought to clarify whether Gb’Socs2 or Gb’Socs5 could suppress cell proliferation. The relative ratio of Gb’cyclinE mRNA in the RLs at 5 dpa was increased in Gb’Socs2^RNAi nymphs (179±20%, n=15, P<0.01; Fig. 3B) compared with control nymphs (100±12%, n=15) as assessed by q-PCR, suggesting that Gb’Socs2 suppresses cell proliferation in the regenerating legs, perhaps by repressing the overactivation of JAK/STAT signalling.

In conclusion, we provide evidence that JAK/STAT signalling promotes blastema formation and that SOCS negatively regulates JAK/STAT signalling activity during tissue regeneration. Taken together, these results indicate that our comparative transcriptome analysis in the cricket system is effective in identifying genes that are functionally involved in tissue regeneration.

We thank Itsuro Sugimura at Hokkaido System Science Co., Ltd and Yasuko Kadomura for their technical support.