Selection and Evaluation of Tissue Specific Reference Genes in Lucilia sericata during an Immune Challenge

The larvae of the common green bottle fly Lucilia sericata (Diptera: Calliphoridae) have been used for centuries to promote wound healing, but the molecular basis of their antimicrobial, debridement and healing functions remains largely unknown. The analysis of differential gene expression in specific larval tissues before and after immune challenge could be used to identify key molecular factors, but the most sensitive and reproducible method qRT-PCR requires validated reference genes. We therefore selected 10 candidate reference genes encoding products from different functional classes (18S rRNA, 28S rRNA, actin, β-tubulin, RPS3, RPLP0, EF1α, PKA, GAPDH and GST1). Two widely applied algorithms (GeNorm and Normfinder) were used to analyze reference gene candidates in different larval tissues associated with secretion, digestion, and antimicrobial activity (midgut, hindgut, salivary glands, crop and fat body). The Gram-negative bacterium Pseudomonas aeruginosa was then used to boost the larval immune system and the stability of reference gene expression was tested in comparison to three immune genes (lucimycin, defensin-1 and attacin-2), which target different pathogen classes. We observed no differential expression of the antifungal peptide lucimycin, whereas the representative targeting Gram-positive bacteria (defensin-1) was upregulated in salivary glands, crop, nerve ganglion and reached its maximum in fat body (up to 300-fold). The strongest upregulation in all immune challenged tissues (over 50,000-fold induction in the fat body) was monitored for attacin-2, the representative targeting Gram-negative bacteria. Here we identified and validated a set of reference genes that allows the accurate normalization of gene expression in specific tissues of L. sericata after immune challenge.


Introduction
Maggot debridement therapy (MDT) was established at the beginning of the 20 th century and became a popular treatment for chronic and recalcitrant wounds. The advent of antibiotics made MDT largely obsolete [1], but intractable wounds still present a challenge for modern medicine and MDT has re-emerged as an alternative therapeutic approach [2,3]. One hallmark of the renewed popularity of MDT is its approval as a medical device by the United States Food and Drug Administration in 2004 (FDA case number K033391), and sterile larvae of the common green bottle fly Lucilia sericata (Diptera: Calliphoridae) have therefore become the exclusive species used for MDT.
Sterile maggots are applied to a wound bed where they debride the necrotic tissue, disinfect the wound and promote wound healing [4,5]. Debridement is accomplished by the secretion of enzymes such as collagenases [6], proteases [7] and nucleases [8], which liquefy the tissue and allow the uptake of nutrients by the maggots. Some of these enzymes can also inhibit or degrade biofilms [8,9,10].
Disinfection is facilitated by the uptake and digestion of microorganisms [11], and by the secretion of antimicrobial metabolites [12,13,14] (e.g. Seraticin) and a variety of antimicrobial peptides (AMPs) [15,16,17,18]. These molecules have the potential to address the emergence and spread of antibiotic resistant pathogens [19]. Indeed, maggot secretions have shown promising activities against clinically relevant strains of drug-resistant pathogens and provide a source of lead structures for new antibiotics [13,14,17]. The antimicrobial activity of maggot secretions has a dose-dependent relationship according to the number and nature of bacteria the larvae encounter [20,21,22]. The secretions of sterile larvae do not possess antimicrobial activity whereas larvae confronted with Staphylococcus aureus or Pseudomonas aeruginosa (the most prominent Gram-positive and Gram-negative bacteria in chronic wounds, respectively [23,24,25,26,27,28]) produce distinct antimicrobial secretions that are not mutually antagonistic [20]. Both bacteria are well known to form biofilms, which protect them against external cues and considerably hamper the treatment [29,30,31].
MDT can also stimulate cell growth and wound healing but this is poorly understood. Maggot secretions can suppress pro-inflammatory responses [32,33], induce blood clotting [34], encourage fibroblasts to spread through the wound [35,36] and regulate the activation of human complement system [37]. Medical maggots also secrete allantoin and urea, which are thought to promote wound healing and are included in a number of medical products [38,39].
Although the molecular basis of MDT has been investigated, little is known about the key molecules responsible for these processes. Next-generation sequencing provides a rapid approach for the generation of genome and transcriptome datasets from non-model organisms such as L. sericata. These datasets allow the identification of medically and industrially relevant genes. To precisely characterize these genes based on their differing expression profiles in specific tissues or under specific conditions other approaches are necessary. The quantitative reverse transcription real-time polymerase chain reaction (qRT-PCR) is a powerful and sensitive method for the quantitation of gene expression [40,41] but is highly dependent on appropriate normalization to correct systemic variations [42,43,44]. The most reliable normalization is achieved by comparing gene expression profiles to stably-expressed reference genes [43,45,46]. A universal reference gene that remains stable under all conditions is unlikely to exist, so panels of candidate reference genes must be validated for certain experimental settings, such as stability during an immune challenge experiment. Here we investigated a panel of 10 candidate L. sericata reference genes in different tissues to allow the normalization of gene expression data in immune challenge experiments, to validate the identification of genes involved in maggot therapy.
L. sericata larvae were obtained from BioMonde GmbH (Barsbüttel, Germany). First instar larvae were reared on Columbia Agar with Sheep Blood PLUS (Thermo Scientific Oxoid) for 72 hours (h) at 28°C in the dark. The fed larvae were cleaned using sterile water before immune challenge or dissection.

Immune challenge and zone of inhibition assays
Larvae were placed in Petri dishes on ice to reduce motility and facilitate injection of larvae with Pseudomonas aeruginosa. The dorsal posterior was pricked with a sterile needle (wounded) or with a needle dipped in P. aeruginosa (DSM 50071) suspension (OD 600 = 60) in phosphate buffered saline (immune-challenged). Treated larvae were supplied with fresh blood agar and incubated for additional 24 h at 28°C in the dark followed by the dissection of individual tissues. For the zone of inhibition assays, 7 ml of 1% LB agar per plate was cooled to 42°C, supplemented with 7 μl of fresh E. coli D31 [47] culture (OD 600 = 0.5) and poured into a Petri dish, before 3-mm wells were stamped into the agar using a sterile hole puncher. After 24 h equal volumes (3 μl) of hemolymph from naïve, wounded and immune-challenged larvae were collected and immediately applied to the agar wells. The plates were incubated at 37°C for 24 h and 3 μl of 100 mg/ml ampicillin were used as a positive control.

Tissue dissection, RNA isolation and cDNA synthesis
Larvae were cooled on ice and dissected by ripping the dorso-anterior cuticle in sterile DEPCtreated PBS under a binocular microscope. The midgut, hindgut, salivary glands, crop, fat body and nerve ganglion were harvested following Freeman and Bracegirdle [48] and equivalent amounts of tissue-specific RNA were collected by preparing three pools of midgut (n = 5), hindgut (n = 10), salivary glands, crop, fat body and nerve ganglion (each n = 25) in RA1 buffer (Macherey-Nagel) before RNA isolation.
Total RNA from naïve and immune-challenged larvae (n = 5) subsequently pooled in equimolar RNA amounts, as well as the dissected larval tissues listed above, were extracted using the NucleoSpin RNA kit (Macherey-Nagel) including a 15-min DNase I on-column digest. The concentration, purity and quality of the RNA were determined by spectrophotometry (Take3, BioTek) and agarose gel electrophoresis (S1 Fig). Only samples with A 260 /A 280 and A 260 /A 230 > 1.8 and at least one sharp band representing 18S rRNA were used for cDNA synthesis. RNA samples that did not meet these criteria were cleaned and concentrated by sodium acetate precipitation [49].
Complementary DNA was synthesized using 1.5 μg of total RNA, oligo(dT) 18 primers and the First Strand cDNA Synthesis Kit (Thermo Fisher Scientific) according to manufacturer's recommendations. The resulting cDNA was diluted to a working concentration of 400 pg/μl, divided into aliquots and stored at -80°C.
Candidate gene sequence assembly L. sericata reference gene candidates were selected based on known reference genes from the closely-related species Lucilia cuprina [50] and other arthropods [51,52]. The peptide sequences of each gene (or the nucleotide sequences of 18S and 28S rRNA) were queried against the L. sericata transcriptome database [15] using the BLAST algorithm. BLASTn was used to find all reads matching the search results and the pool of these reads was expanded to include further paired-end read partners. The final group of reads was assembled using Trinity [53] and this step was repeated until the complete coding sequences were acquired. Finally, Gap5 [54] was used to verify the finished sequences with mapped reads assembled by Bowtie [55].

Primer design and evaluation
Gene-specific primers were designed using Oligo Explorer v1.1.2 (http://www.softpedia.com/ get/Science-CAD/Oligo-Explorer.shtml) to yield primers 19-23 nucleotides in length with amplification products of 50-210 bp and T m values of~60°C. All primers were subsequently tested in a standard curve assay including melt curve analysis using the StepOnePlus Real-Time PCR System (Applied Biosystems). The cDNA concentration ranging from 50 ng to 3.2 pg of total RNA was used in 5-fold dilutions with the cycling conditions recommended by the manufacturer. Briefly, hot-start PCR with denaturation at 95°C for 10 min was followed by 40 cycles of 95°C for 15 s and 60°C for 60 s, and finally the melt curve analysis with a temperature increase from 60°C to 95°C in 0.5°C steps. Reaction efficiency was calculated using StepOne Software v2.3 and only primers with an efficiency of 90-110%, R 2 0.99 and a single sharp melt curve peak were used for reference gene evaluation.

Quantitative real-time RT-PCR
Quantitative real-time RT-PCR was carried out on a StepOnePlus system (Applied Biosystems) using optical 96-well plates (Applied Biosystems). The total reaction volume of 10 μl contained 5 μl Power SYBR Green PCR Master Mix (Applied Biosystems), 1 μl (400 pg) cDNA and 300 nM of each primer except for β-tubulin (150 nM), EF1α (150 nM) and actin (450 nM forward primer and 150 nM reverse primer). All reactions were carried out in three technical replicates under the reaction conditions stated above. Baseline correction was performed automatically by StepOne Software v2.3 and the quantification cycle (Cq) was always determined at an intensity of 0.15. All cDNA samples and corresponding RNA without reverse transcription (no-RT control) were tested with the 18S rRNA primers to estimate the remnants of genomic DNA [56] and only samples with a ΔCq 10 were accepted.

Normfinder and GeNorm analysis
The Normfinder and GeNorm algorithms were used for reference gene assessment. Normfinder [57] was used as an add-in for Microsoft Excel (http://moma.dk/normfinder-software), and GeNorm [46] was used as part of the R package NormqPCR with R version 3.1.1 [58]. In both cases, raw Cq values were used and log transfer was performed by the software.

Analysis of gene expression following immune challenge
Gene-specific primers were designed for the coding sequence of three L. sericata immune genes lucimycin, defensin-1 and attacin-2, which are differentially expressed when the larval immune system is challenged [15]. Genes were chosen based on their activity spectrum, with lucimycin targeting fungi [16], defensin-1 targeting Gram-positive [59] and attacin-2 targeting Gram-negative bacteria [59], respectively. Gene expression profiles upon immune challenge were determined by qRT-PCR using the ΔΔCq method [60] in Rest2009 (http://www.genequantification.de/rest-2009.html). The data were normalized using three different combinations of genes from the reference gene assessment: (a) the two best reference genes for every sample (GeNorm); (b) the three overall best reference genes (GeNorm and Normfinder); and (c) the six overall best reference genes (GeNorm). Efficiencies greater than 100% were set to 100% because Rest2009 does not allow higher values.

Selection of candidate reference genes
We monitored the expression level of 10 L. sericata candidate reference genes before and after challenging the larvae with P. aeruginosa (Table 1). Our goal was to identify reference genes appropriate for different larval tissues, with the key criterion that gene expression is not affected by immune challenge. The tissues were selected based on their role in secretion and/or digestion (salivary glands, crop, midgut and hindgut) and supplemented by additional tissues (nerve ganglion and fat body). Candidates were chosen based on earlier studies in arthropods using orthologs of these genes [50,52]. To minimize the risk of co-regulation, we selected genes representing different functional classes of proteins. The L. sericata transcriptome database was explored to identify orthologs and the coding sequences of the 10 selected candidate reference genes were submitted to GenBank (Table 1).

Quantitative RT-PCR
The efficiencies of each qRT-PCR primer pair were generally high. Based on the standard curve slopes determined using StepOne software, the efficiencies ranged from 90% to 109% with R 2 values > 0.99 for all pairs (Table 1). Linear behavior was observed over a concentration range spanning four orders of magnitude (50 ng to 3.2 pg of cDNA). The absence of primer dimers was confirmed by a melting curve analysis, which resulted in only one sharp peak for each amplicon under our experimental conditions (

Zone of inhibition assays and initial analysis of attacin-2 expression
To monitor the changes in larval reference gene expression triggered by interaction with bacteria, we chose the opportunistic pathogen P. aeruginosa and direct infection by pricking the L. sericata larvae with a bacteria-coated needle [61]. P. aeruginosa has previously been used to induce an immune response [15] and is also a common pathogen found in chronic wounds. Two independent tests were used to monitor the success of the immune challenge.
The first was a zone of inhibition assay with E. coli strain D31 [47]. Hemolymph from naïve, wounded and immune-challenged larvae were applied to the agar wells and the plates were incubated at 37°C for 24 h. Only the hemolymph from immune-challenged larvae produced distinct inhibition zones (Fig 2a). As a second test, we carried out an initial analysis of gene expression profiles to monitor changes at the RNA level triggered by immune challenge, using attacin-2 mRNA as a marker because this AMP is known to be upregulated under similar conditions [15]. The raw Cq values for attacin-2 mRNA from two distinct populations, i.e. the naïve samples (Cq~30) and the immune-challenged (Cq~20), are shown in (Fig 2b). No naïve sample showed a lower value than the corresponding immune-challenged sample. These raw data supported our initial hypothesis that L. sericata attacin-2 mRNA is induced by immune challenge with a Gram-negative bacterium, but can only be validated by normalizing against the expression levels of reliable reference gene(s).

Validation of candidate reference genes using Normfinder
The model-based Normfinder approach provides a direct estimation of variance in so-called stability values, with lower numbers equating a more stable expression. The program does not sequentially eliminate unsuitable genes so the pre-exclusion of genes with large variances is necessary to achieve reliable results.
Normfinder analysis of the 10 candidate reference genes (Table 2) ranked RPS3, RPLP0 and EF1α as the three most stable genes considering immune challenge. The worst stability values (approximately 2-fold higher across all samples) were observed for 18S rRNA, actin and 28S rRNA. The overall ranking was not transferable to all specific tissues or to whole larvae because the stability of the candidates differed according to the sample type. Intragroup variation was generally low with only 28S rRNA showing variation (S1 Table). Intergroup variation was generally much higher, with the highest scores generated by 18S rRNA, 28S rRNA and actin (S2 Table).

Validation of candidate reference genes using GeNorm
The GeNorm algorithm calculates the pairwise variation among all tested genes and assigns stability measures (M). In every cycle, the gene with the highest stability measure (i.e. the least stable) is excluded until only the two best genes remain. GeNorm analysis of our 10 candidate reference genes (Table 3) yielded similar results to Normfinder. GeNorm ranked RPS3, RPLP0 and EF1α as the three best genes and 18S rRNA, actin and 28S rRNA as the three worst, with M values of 2 to 4-fold higher than the best candidates across all samples. As for the Normfinder data, the overall GeNorm ranking list was not transferable to all individual tissues or whole larvae.
To illustrate the differences in stability measures we also included attacin-2, which is strongly induced by immune-challenge and therefore expected to be unstable. As anticipated, the M value for attacin-2 was~20-fold higher than the most stable candidate genes in each sample (Table 3).
Because neither algorithm selected a single reference gene suitable for all samples, we used GeNorm pairwise variation comparison and select the most appropriate number of reference genes for normalization, which is an appropriate strategy when single reference genes are unsuitable [46]. This comparison (V n/n+1 ) examines the normalization factors NF n and NF n+1 in each analysis cycle. If the calculated value falls below a set threshold of 0.15 [46], the addition of further reference genes does not improve the quality of normalization. As shown in Fig 3, we found that most tissues reached the recommended threshold of 0.15 with just two reference genes, whereas three reference genes were required for the normalization of expression levels in the salivary glands (V 2/3 = 0.17). When comparing overall variation, the 0.15 threshold is reached at V 6/7 which means that six reference genes would be necessary for appropriate simultaneous normalization within all tested samples.
The Cq values for lucimycin, defensin-1 and attacin-2 mRNA were therefore normalized using three combinations of candidate reference genes: (a) the two best candidate reference genes for every sample defined by GeNorm analysis; (b) the overall best three reference genes defined by both algorithms (RPS3, RPLP0 and EF1α); and (c) the overall best six reference genes defined by GeNorm pairwise variation comparison. Our results show that the antifungal peptide lucimycin is not differentially expressed (Fig 4a), whereas defensin-1, which is specifically targeting Gram-positive bacteria, shows different upregulation in salivary glands, crop, nerve ganglion and fat body, ranging from a 8-fold induction in the nerve ganglion to a 290-fold induction in the fat body (Fig 4b). The strongest upregulation was monitored for attacin-2 that targets primarily Gram-negative bacteria. Attacin-2 was strongly induced by immune challenge in all tested tissues, ranging from a 37-fold induction in the salivary glands to a 51,463-fold induction in the fat body (Fig 4c). As shown in Fig 4, all three reference gene combinations lead to similar expression results with slight variations in relative values.

Discussion
MDT has undergone a renaissance since the end of the 20 th century, but although recent studies have addressed the underlying mechanisms many of the key proteins responsible for the beneficial effects of medical maggots remain unknown [4,5]. The falling cost of next generation sequencing means that the sequencing of genome and transcriptome datasets from non-model organisms can now be used as a first step to screen for medically and industrially relevant sequences. In a second step qRT-PCR can validate these candidates, but only if appropriate reference genes are available for the normalization of expression data [43,45,46]. This approach is ideal for the identification of L. sericata genes with important roles in wound healing. Thus such genes are likely to show different expression profiles in specific tissues or are modulated in the presence of wound pathogens. We investigated the tissue specific expression stability of 10 candidate reference genes and 3 immune genes (lucimycin, defensin-1 and attacin-2), which target different pathogen classes. We monitored the expression of these immune genes in different L. sericata tissues before and after the immune challenge with the Gram-negative bacterium P. aeruginosa, which was introduced into third-instar L. sericata larvae by pricking with a contaminated needle. We chose this invasive strategy to ensure that the comparator underwent a strong change in expression allowing the stability of the candidate reference genes to be tested robustly, even though this far exceeds the intensity of the bacterial challenge that maggots would encounter in human wounds. As expected, we monitor the strong expression of attacin-2 (targeting Gram-negative bacteria [59]) in all larval tissues reaching its maximum in fat body (over 50,000-fold). We also did not observe any differential expression of the antifungal peptide lucimycin, which is known to have no effect against Gram-negative bacteria [16].
Interestingly, defensin-1 as a representative targeting Gram-positive bacteria [59], was differentially upregulated in salivary glands, crop, nerve ganglion and fat body. This unexpected result can be easily explained by the wounding procedure itself, which has been shown to induce insect immune response previously [62]. Our results demonstrate that even under such harsh conditions, which heavily perturbed the system and led to the strong induction of  attacin-2, the expression of most of the candidate reference genes remained stable. This success rate reflects our choice of genes that have been used as references in other insects, allowing the knowledge-based pre-selection of candidates that are likely to be suitable as reference genes in L. sericata [50,63,64,65,66,67]. The Normfinder and GeNorm algorithms were unable to identify a single optimal reference gene or an ideal combination of two to three reference genes that worked consistently across tested samples. RPLP0, EF1α and RPS3 were identified as the most stable reference genes across all samples even though each of them ranked low for stability in at least one sample. They are all involved in protein expression [68,69,70] and co-regulation cannot be ruled out as explanation for their high ranking. However, comparison with the other candidate genes does not indicate generally co-regulated expression. Although it may be necessary to use specific candidates for particular tissues to achieve the most accurate data normalization, the use of these overall three best reference genes should be more than adequate for most experimental scenarios, as illustrated by our expression analysis of three immune genes (Fig 4).
Among the remaining candidates, some (e.g. GST1 and GAPDH) occupy mid-ranking positions in many tissues and have poor scores in others but never rank as the most suitable reference gene, whereas others (e.g. PKA and β-tubulin) also occupy mid-raking position in many tissues but are the most stable reference genes in others, i.e. the salivary glands and the nerve ganglion (PKA) or the crop (β-tubulin). Their overall mid-to-low ranking reflects their tendency to show good scores in a few tissues but poor scores in most others.
Remarkably, actin has a low overall score but good scores in several tissues (second most stable reference in three tissues and third in another). However, the standard deviation of 2.1 Cq is the highest among all the tested candidate genes, reflecting high intergroup variation (S1 Table) rather than changes in expression induced by the immune challenge, which would lead to higher intragroup variations. It is not surprising that a fundamental structural protein like actin would be expressed at different levels in diverse tissues such as the gut and the fat body. Therefore, despite its low overall ranking, actin would be useful as a reference gene in most individual tissues but not in more diverse collections of samples.
Compared to the other candidates, the genes for 18S and 28S rRNA were the least stable and therefore the least suitable for normalization, which contrasts with the findings in other arthropods [50,71]. This may be a species-dependent phenomenon, but may also reflect the technical approach, e.g. the method used for RNA isolation and cDNA synthesis. For example, we used oligo(dT) 18 primers which only partly reverse transcribe the rRNA genes due to the lack of a canonical polyadenylate tail.
While dissecting the fat body from immune-challenged larvae we observed the remodeling and/or degradation of this tissue. The physical appearance of the fat body changed from clusters of round cells to loosely-associated amorphous cells which yielded lower amounts of total RNA with poorer quality compared to the other tissues. After immune challenge, fat body total RNA from all biological replicates had to be precipitated and concentrated to meet the required criteria for cDNA synthesis. The insect fat body is responsible for AMP production [72] and two defense pathways have been described in Drosophila melanogaster, i.e. the Toll pathway against Gram-positive species and the immune deficiency (IMD) pathway against Gram-negative species [73]. Constitutive triggering of the IMD pathway has been shown to induce apoptosis [74]. Our results show that although attacin-2 is expressed ubiquitously and globally upregulated by immune challenge, the strongest induction (over 50,000-fold) occurs in the fat body. We used P. aeruginosa, which as a Gram-negative species is targeted by attacins [59] and we hypothesize that the IMD pathway plays a major role in this immune response. The induction of attacin-2 may therefore be part of this IMD pathway response, ultimately resulting in apoptosis which may explain the physical changes in the fat body we observed.
There is probably no universally ideal reference gene in any experimental system so it is important to distinguish between the best gene and combination of genes suitable for different experimental settings. Despite the diverse expression profiles of the 10 candidate reference genes in our dissected tissue samples, several combinations of these genes gave identical results when used for the normalization of immune genes expression during an immune challenge, suggesting that the analysis of six reference genes per sample or different reference genes in each tissue is unnecessary. From our own experience, the reliable normalization of L. sericata expression data among all dissected tissues (Fig 4) is also achieved using the three best reference genes based on our overall rankings (RPLP0, EF1α and RPS3), which balances accuracy with convenience and the cost of additional experiments. The primer pairs and reference gene candidates evaluated in this study provided a valuable tool for the normalization of gene expression data in medical maggots, which will facilitate the identification and functional analysis of genes which are responsible for the beneficial effects of maggot therapy.