Comparative Transcriptome Analysis of Bombyx mori (Lepidoptera) Larval Midgut Response to BmNPV in Susceptible and Near-Isogenic Resistant Strains

Bombyx mori nucleopolyhedrovirus (BmNPV) is one of the primary pathogens causing severe economic losses in sericulture. However, the molecular mechanism of silkworm resistance to BmNPV remains largely unknown. Here, the recurrent parent P50 (susceptible strain) and the near-isogenic line BC9 (resistance strain) were used in a comparative transcriptome study examining the response to infection with BmNPV. A total of 14,300 unigenes were obtained from two different resistant strains; of these, 869 differentially expressed genes (DEGs) were identified after comparing the four transcriptomes. Many DEGs associated with protein metabolism, cytoskeleton, and apoptosis may be involved in the host response to BmNPV infection. Moreover, some immunity related genes were also altered following BmNPV infection. Specifically, after removing genetic background and individual immune stress response genes, 22 genes were found to be potentially involved in repressing BmNPV infection. These genes were related to transport, virus replication, intracellular innate immune, and apoptosis. Our study provided an overview of the molecular mechanism of silkworm resistance to BmNPV infection and laid a foundation for controlling BmNPV in the future.


Introduction
The silkworm, Bombyx mori L. (Lepidoptera: Bombycidae) has been domesticated for production of cocoons for more than 5000 years. Silkworm rearing and the silk industry still play an important role in China, India and many other developing countries. B. mori is also a good model for the study of insect genetics and immunology [1][2][3][4]. Bombyx mori nucleopolyhedrovirus (BmNPV) is the principal silkworm pathogen and causes serious economic losses in sericulture every year. Among numerous silkworm strains, most are susceptible to BmNPV infection, although a few strains exhibit resistance [5]. The heredity of silkworm resistance against BmNPV infection is a relatively complicated process because resistance is controlled both by major dominant genes and multiple genes of micro-effect [6].
A series of studies have made significant progress in understanding silkworm resistance against BmNPV infection. Xu et al. reported that Bms3a was potentially involved in resistance to BmNPV infection [7,8]. B. mori lipase-1, serine protease-2 and alkaline trypsin protein extracted from the digestive juice of larvae midguts showed strong antiviral activity in vitro [9][10][11]. Using comparative proteomics, arginine kinase was found to be involved in the antiviral process of different resistant strains of silkworm [12]. In our laboratory, a total of 12 proteins that are potentially involved in viral infection were identified using one-and two-dimensional electrophoresis followed by virus overlay assays. These proteins could be categorized into the following groups: endocytosis, intracellular transportation, and host responses [13]. Immune responses were found to be synergistically regulated by the Toll, Janus kinase/signal transducers and activators of the transcription (JAK/STAT) and immune deficiency (IMD) pathways, which could act as an important defense against exogenous pathogenic infection in conjunction with subsistent pathogen recognition receptors and response proteins [14][15][16][17][18]. However, the molecular mechanisms of silkworm resistance to BmNPV infection are still not fully elucidated.
In recent years, the high-throughput nature of next generation sequencing (NGS), using platforms such as Illumina HiSeq TM 2500 have provided fascinating opportunities in the life sciences and dramatically improved the efficiency and speed of gene discovery, especially in the research of host cell responses to exogenous pathogenic infection [19]. For example, Hu et al. obtained numerous differentially expressed genes (DEGs) involved in metabolism, immunity, and inflammatory responses in Microtus fortis following infection with Schistosoma japonicum based on comparative transcriptome analysis [20]. Diege et al. examined different fish tissues infected with salmon anemia virus (ISAV) using high-throughput transcriptomics and found a strong correlation between functional modules and viral-segment transcription [16]. NGS technology was also used to explore the molecular mechanism of silkworm resistance against exogenous pathogens. Kolliopoulou et al. reported that several genes related to physical barriers, immune response, proteolytic/metabolic enzymes, heat-shock proteins, and hormonal signaling were possibly involved in silkworm resistance against Bombyx mori cytoplasmic polyhedrosis virus (BmCPV) infection; although these genes might be induced by the virus in order to increase infectivity [21]. Additionally, several candidate genes, such as BmEts, BmToll10-3 and Hsp20-1, have been identified in the initial stage of BmNPV infection by analyzing the global transcriptional profile of silkworm cell lines and heads following BmNPV infection [22,23].
In order to gain a global view of the molecular changes in silkworms during BmNPV infection, we selected near-isogenic line BC9 and recurrent parent P50 for transcriptome sequencing. Through comparative analysis of the transcriptomes from these two strains, a total of 869 DEGs were obtained, which included many genes potentially related to BmNPV-resistance. Our results may provide some reliable evidence to clarify the BmNPV-resistance molecular mechanism in silkworm.

Virus and Silkworm
BmNPV (T3 strain) was maintained in the Key Laboratory of Sericulture, Anhui Agricultural University, Hefei, China. The virus was obtained from the hemolymph of infected larvae and purified by repeated and differential centrifugation according to the protocol developed by Rahman et al. [24]. The concentration of the virus (OB/mL) was determined by hemocytometer.
The recurrent parent P50 (susceptible strain), the donor parent A35 (resistant strain) and the near-isogenic line BC9 were maintained in our laboratory. The near-isogenic line was constructed according to the protocol used by Yao et al. [6]. In brief, the recurrent parent P50 was crossed to the donor parent A35; progeny were repeatedly backcrossed with the recurrent parent for nine generations and each progeny was screened by BmNPV.
The first three instar larvae were reared on a fresh artificial diet at 26±1°C, 75±5% relative humidity, and a 12 hours day/night cycle. The rearing temperature for the last two instars was reduced to 24±1°C; other conditions were unchanged. On the first day of fifth instar, all larvae were starved for 24 hours and then fed with 5 μL BmNPV suspended in sterile water (1.0×10 5 OB/mL) per larva orally; the control group was treated with sterile water. BmNPV occlusion bodies (OB) began fast proliferation at approximately 24 hours post inoculation (hpi) [25]; thus, this time was considered optimal for sample collection. Silkworm larvae were dissected and the midgut tissues were removed and then washed in PBS (137 mM NaCl, 2.7 mM KCl, 4.3 mM Na 2 HPO 4 , and 1.4 mM KH 2 PO 4 , pH 7.4) prepared with diethy pyrocarbonate (DEPC) (Sangon, China) treated H 2 O. Thirty larvae midguts were mixed together to minimize individual genetic differences. Samples were flash-frozen in liquid nitrogen and pulverized, and 100 mg of sample were added directly into a RNAase free microcentrifuge tube containing 1.0 mL TRIzol Reagent (Invitrogen, USA) and stored at -80°C for later use.

Silkworm strain resistance level bioassays
The level of silkworm resistance to BmNPV was tested following the protocol developed by Cheng et al. [26]. The fourth instar larvae were inoculated with BmNPV at different concentrations; inoculations were conducted in triplicate. The level of silkworm resistance was calculated using IBM SPSS Statistics 20 (IBM, USA).

RNA extraction
The silkworm midguts dissolved in TRIzol Reagent were homogenized. Total RNA of midguts were extracted according to the manufacturer's protocol. Concentrations were quantified using a NanoDrop 2000 spectrophotometer (Thermo Scientific, USA). The purity of all RNA samples were assessed at an absorbance ratio of A 260/280 and A 260/230 , and the integrity of the RNA was confirmed by 1% agarose gel electrophoresis.
Library construction, Illumina sequencing and read assembly Fragment interruption, cDNA synthesis, addition of adapters, PCR amplification and RNA--Seq were performed by Beijing BioMarker Technologies (Beijing, China). The standard Illumina methods and protocols were adopted to prepare and sequence the cDNA libraries. NEBNext Poly(A) mRNA Magnetic Isolation Module (NEB, USA) was used to enrich mRNA, and then the cDNA library was constructed using the NEBNext mRNA Library Prep Master Mix Set for Illumina (NEB, USA) and NEBNext Multiplex Oligos for Illumina (NEB, USA). The size of the library insert fragments was determined by 1.8% agarose gel electrophoresis, and the fragments were quantified using a Library Quantification Kit/Illumina GA Universal (Kapa, USA). Suitable fragments were selected as templates and sequenced on an Illumina HiSeq TM 2500 using paired-end technology. Three biological replicates were used to minimize sample differences.
In order to obtain clean and high-quality reads for sequence assembly, the raw reads were filtered by removing adaptor sequences, low-quality sequences (reads with ambiguous bases 'N') and reads with 10% > Q < 20% bases [27,28]. The Trinity assemble program was used to assemble the clean reads into contigs, which covered more full-length transcripts through a broad range of expression levels [29]. The resultant contigs were added to transcripts based on paired-end information. The longest transcript from alternative splicing transcripts was selected as the unigene. These unigenes were combined to produce the final assembly and used for annotation.

Functional annotation
To annotate unigenes, different sequences were searched by BLASTx against the NCBI nonredundant protein (nr) database and other databases, including Swiss-Prot protein database, the Kyoto Encyclopedia of Genes and Genomes (KEGG) and the Cluster of Orthologous Groups (COG) database. Gene Ontology (GO) annotations, including molecular functions, biological processes, and cellular components, were obtained using the Blast2GO program (https://www.blast2go.com/) [30,31]. All searches were performed with an E-value < 10 −5 . Fragments per kilobase of transcript per million fragments mapped (FPKM) was calculated to represent the expression abundance of the unigenes [32]. FPKM may reflect the molar concentration of a transcript by normalizing for RNA length and for the total read number.

Differentially expressed genes (DEGs) analysis
After normalizing genes expression levels, DEGs were obtained by pair-wise comparison of the four transcriptome libraries using IDEG6 software [33]. An FPKM fold change of 1.2 or 0.83 between two libraries was defined as the reference standard, with the Benjamini-Hochberg false discovery rate (FDR < 0.01) used to adjust the p-values.

Real-time quantitative PCR analysis
In order to validate the results from our transcriptome sequencing analysis, the relative expression levels of 15 randomly selected genes were confirmed by reverse transcription quantitative PCR (RT-qPCR). Additionally, 9 genes with well reported previously were selected to further validate the genes of interest that might be involved in BmNPV resistance. All the Primers are listed in Table 1. RT-qPCR reactions were prepared with the SYBR Premix Ex Taq TM Kit (TaKaRa), following the manufacturer's instruction. Reactions were carried out in Bio-Rad CFX96TM Real-Time System (Bio-Rad, USA). The thermal cycling profile consisted of an initial denaturation at 95°C for 30 s, 40 cycles at 95°C for 5 s, and 60°C for 30 s. All samples were performed in triplicate. Relative expression levels were calculated using the 2 -44Ct method following the protocol of Livak et al. [34]. In this study, B. mori ribosomal protein s3 (BmRPS3) gene was used as a reference gene. Statistical analysis was conducted using the SPSS software (IBM, www.ibm.com).

BmNPV infectivity in different silkworm strains
The LC 50 value was used to evaluate the resistant level of silkworm to BmNPV infection. The LC 50 value of A35 was approximately 26-fold greater than that of BC9 and over 500-fold greater than that of P50. The value of BC9 was 23-fold greater than that of P50 ( Table 2).

Overview of the silkworm transcriptome
Transcriptome sequencing is an efficient technology for comparing gene expression levels in different samples, and in our study, it was used to search and analyze DEGs among P50, BC9 control and treatment groups. A total of four cDNA libraries were sequenced: P50-(P50 treated with sterile water), P50+ (P50 infected with BmNPV), BC9-(BC9 treated with sterile water), BC9+ (BC9 infected with BmNPV), with each group created in triplicate. After removing the adaptors and low quality sequences, 144,439,382 sequence reads were obtained ( Table 3). The GC content of each of the four libraries was approximately 50%, and CycleQ30% was greater than 89.91% for each library. Thus, the quality and accuracy of the sequencing data were sufficient for further analysis. Most of the reads matched silkworm genomic locations. All the unigenes matched previously described sequences with approximately 70% coverage. The length distribution of unigenes had similar patterns among the four libraries, suggesting that there was little bias in the construction of the four cDNA libraries (Fig 1).

Unigenes annotation and classification
In order to annotate the unigenes, reference sequences were searched using BLASTX against the NCBI nr database (E-value < 10 −5 ). A total of 12,591 of 13,342 unigenes provided a BLAST result (S1 Table). S2 Table shows the species with the closest match for each unigene. Most of Table 1. Primers used in RT-qPCR for validation of DEGs.

RT-qPCR validation of differentially expressed transcripts
In order to determine the reliability of the transcriptome sequencing, the relative expression levels of 15 randomly selected genes were analyzed by RT-qPCR (Fig 2A). The results were consistent with the transcriptome data. For example, the gene encoding peptidoglycan-recognition protein was down-regulated in both RNA-seq and RT-qPCR analyses, with a similar fold change. The lipase 1 gene was significantly up-regulated in the resistant strain BC9 after  infection with BmNPV, which was also consistent with previously reported results [10]. Linear regression analysis of the correlation between RT-qPCR and RAN-seq ( Fig 2B) showed an R 2 (goodness of fit) value of 0.9169 and a corresponding slope of 1.5281, suggesting a strong positive correlation between RT-qPCR and transcriptome data. Therefore, the transcriptome data were satisfied for further analysis.

Differentially expressed genes (DEGs) and their possible roles in host response to BmNPV
To further elucidate which DEGs had a potential role in antiviral response, a Venn diagram was constructed. A total of 285 DEGs were found to be differentially regulated when comparing P50+ and P50-, of which 122 were up-regulated and 163 were down-regulated. Similarly, 193 DEGs were found to be differentially regulated in the comparison of BC9+ vs. BC9-, with 56 genes up-regulated and 137 down-regulated. In addition, there were 154 DEGs differentially regulated in the comparison of BC9-and P50-, among which 78 genes were up-regulated and 76 genes were down-regulated (Fig 3, Table 4, S3 Table). There were 197, 119, 82 unique DEGs in P50+ vs. P50-, BC9+ vs. BC9-and BC9-vs. P50-, respectively. GO assignments were used to assign a functional classification to these DEGs. For cellular components, the number of unique DEGs fell into the macromolecular complex classification was distinct in BC9 and P50 following BmNPV infection. For molecular functions, the number of transporter activity related unique DEGs was distinct in BC9 and P50 following BmNPV infection. For biological progresses, the number of unique DEGs involved in metabolism processes and localization was distinct in BC9 and P50 following BmNPV infection (Fig 4).

Analysis of DEGs associated with protein metabolism, cytoskeleton, and apoptosis
The comparisons of BC9+ vs. BC9-and P50+ vs. P50-identified many DEGs that might either be involved in silkworm defense against BmNPV or facilitate BmNPV infection. These genes could be divided into three categories: protein metabolism, cytoskeleton, and apoptosis.
Most of the DEGs (70%) associated with protein metabolism were down-regulated in BC9 following BmNPV infection. In contrast, 80% of the genes in P50 were up-regulated after BmNPV infection. Most of the DEGs (87.5%) associated with the cytoskeleton were up-regulated in P50 after BmNPV infection, but the number of up-regulated genes (50%) decreased in BC9. A total of 19 DEGs associated with apoptosis were identified following BmNPV infection. Ten of these DEGs were altered in BC9 after BmNPV infection. The other DEGs were all overexpressed in P50 following BmNPV infection (Table 5).

Alteration in immunity related gene expression after BmNPV infection in different resistant strains
Pathogen infection stimulated both cellular and humoral responses of insects [33][34][35][36]. Genes participating in innate immunity pathways were identified and analyzed in regards to their potential role in BmNPV infection (Table 6). Thirty DEGs were identified, which could be classified into Toll pathway, Imd pathway, polyphenol oxidase (PPO) pathway, pattern recognition receptor, and antimicrobial peptide. As shown in Table 6, 14 (47%) of these genes were downregulated and 9 (30%) were up-regulated in BC9 after infection with BmNPV, while 17 (57%) were down-regulated and 6 (20%) were up-regulated in P50 after BmNPV infection.

Role of selected DEGs in BmNPV resistance
After removing the genetic background and strain specific immune stress response genes, a total of 22 DEGs were identified as possibly being involved in silkworm resistance to BmNPV (Fig 3, Table 7). For BC9-vs. P50-, 13 genes were up-regulated and the rest were down-regulated. However, the 22 genes were all down-regulated in BC9 following BmNPV infection ( Fig  5). Some genes, including prostatic acid phosphatase, protease inhibitor 6, actin cytoskeleton-    regulatory complex protein PAN1, and EF-hand domain-containing protein, were further analyzed by RT-qPCR (Fig 6). After BmNPV infection, the expression levels of 4 genes were down-regulation in BC9 and A35 (resistant strain) (Fig 6), which was consistent with the transcriptome data. Thus, we deduced that the 22 genes were possibly involved in resistance to BmNPV infection. Gene functions fell into the following categories: transport, virus replication, intracellular innate immunity, and apoptosis.

Discussion
Despite the confirmation of an association between many genes and proteins and resistance to BmNPV, the molecular mechanism of antiviral activities was still unclear. Here, transcriptome sequencing was carried out to identify genes related to BmNPV-resistance in silkworm across the genome. By using the near-isogenic line BC9 (resistant strain) and the recurrent parent P50 (susceptible strain) to study silkworm antiviral mechanisms, some DEGs responding to BmNPV infection were successfully identified after comparing infected groups and controls in the two strains.
Protein metabolism, cytoskeleton, and apoptosis may play an important role in host response to BmNPV infection Cellular and metabolic processes will be dramatically changed after viral infection [35]. In our study, several DEGs that participate in protein metabolism were found to be of interest. For example, solute carrier family 12 is involved in transporting extraordinarily diverse solutes [45], and cystathionine gamma-lyase participates in hydrolysis of cystathionine [46]. Viruses may have to rely on cell proteins to accomplish replication in intercellular regions [47], therefore, the downregulation of the genes involved in protein metabolism could inhibit the replication of BmNPV in the host cell. We speculated that the down-regulation of these genes affected virus replication. Cytoskeleton-dependent intracellular transport is an important strategy for transport of viral particles to different destinations [48]. In this study, some cytoskeleton related genes were found to be of interest, including actin cytoskeleton-regulatory complex protein PAN1 and actin-binding protein. These genes related to actin-coupled endocytosis could promote viral transport [49,50]. We speculated that the down-regulation of the cytoskeleton genes might affect BmNPV transport.
Apoptosis plays a vital role in regulating cell response in Lepidopteran insects during viral infections, where larvae use selective apoptosis and subsequent sloughing of the infected cells in the midgut epithelium to resist virus infection [51,52]. In this study, some genes related to activation of apoptosis were found to be of interest, including cytochrome c, inhibitor of caspase-activated DNase, amyloid precursor protein and B-cell lymphoma protein 2. The activity of caspase-activated DNase was blocked by Hepatitis C virus core at physiological levels, resulting in the inhibition of apoptotic cell death [53]. Amyloid precursor protein is a member of several signaling pathways that are involved in abnormal cell cycles, subsequently leading to apoptosis [54]. B-cell lymphoma protein 2 could bind to BH3 domains of various pro-apoptotic regulators to activate apoptosis [55]. The overexpression of cytochrome c in rabies virus showed a decreased pathogenicity in vitro and in vivo [56]. Based on their role in apoptosis activation, hosts need to increase the expression level of these genes to promote apoptosis when exposed to a virus; this supposition explains the up-regulation of genes involved in apoptosis in the transcriptome following BmNPV infection. Additionally, a significantly higher relative transcriptional level of cytochrome c was detected by RT-qPCR in BC9 and A35 (Fig  6), which indicated the up-regulation of this gene could activate apoptosis to repress further BmNPV infection. However, some genes were down-regulated, including cysteine aspartic acid specific protease 9L, protein kinase A, apoptosis-inducing factor, apoptosis signal-regulating kinase 1, TNF-receptor-associated factor 6, TGF-beta-activated kinase 1, and p90 ribosomal S6 kinase. The inhibition of these genes may strongly impair viral infectivity and virus-induced apoptosis [57][58][59][60][61][62].

Changes to immune gene expression following infection with BmNPV
Although some differential expression of immune genes was observed, this might not be considered biologically important due to low expression levels; low expression levels were also found after BmCPV infection [21]. In this study, most immune related genes were down-regulated, with a few exceptions, such as the 2-fold up-regulation of toll-like protein and lysozyme. The expression of several other genes, including 18 wheeler, macrophage mannose receptor 1, and slit homolog 3 protein, were also up-regulated during the BmNPV infection. Unfortunately, the relationship between immune genes and BmNPV remains unclear and requires further study. We presumed that the low expression levels of immunity related genes may be associated with the disruption of the immune system by BmNPV, similar to the pathogenicity of human immunodeficiency virus (HIV).

Multiple genes have potential roles in repressing BmNPV infection
Based on the Venn diagram, 22 genes of interest were identified, all of which were down-regulated in BC9 following BmNPV infection. These genes were grouped based on their functions, as reported in the literature. These groups, including transport, virus replication, intracellular innate immune, and apoptosis, may play an important role in the process of silkworm resistance to BmNPV. Several genes related to virus transport, including lactase-phlorizin hydrolase (LPH) [37], B (0,+)-type amino acid transporter 1 (BAT1) [38,46], actin cytoskeleton-regulatory complex protein PAN1 (PAN1) [63], and otoferlin [41,42] were identified in this study. The expression level of these genes were all down-regulated in BC9 following BmNPV infection (Fig 5). The resistant strain A35, a donor parent, was used to validate our results. The expression level of PAN1 was chosen for further testing by RT-qPCR. Expression levels of PAN1 were down-regulated in both BC9 and A35 following BmNPV infection (Fig 6). The down-regulation of virus transport-related genes could inhibit the transmembrane and intracellular transport of BmNPV thereby preventing further infection. MFS-type transporter (MFS) [39,40], a transmembrane facilitator, was induced to increase intracellular monovalent ion concentrations, which led to lysis and cell death after HIV-1 infection. The down-regulation of this gene in silkworm could potentially block BmNPV infection.
Protease inhibitor 6 is a member of the protein superfamily that contains TIL functional domain. Zhao et al. used genome sequences to demonstrate that the expression level of the TIL superfamily were down-regulated following BmNPV infection [77], which was consistent with our results (Fig 6). Furthermore, peptidoglycan-recognition protein 2, hemolin, facilitated trehalose transporter Tret1, integument esterase 2, defense protein precursor, and antimicrobial protein 6Tox precursor also showed differential expression following BmNPV infection. The relationship of these genes to BmNPV is an important area for further study.
Four DEGs were down-regulated in BC9 following BmNPV infection, including rasresponsive element-binding protein 1, GTPase-activating protein, trypsin alkaline A and peroxidase (Fig 5). Previous studies revealed that these proteins play a role in virus infections. Mdv1-miR-M4 encoded by Marek's disease virus efficiently targeted the 3' untranslated regions of ras-responsive element-binding protein 1 (RREB1) [78]. TBC domain proteins belonging to a GTPase-activating protein were knocked out by double stranded RNA interference (RNAi), which led to a decrease in the level of transcripts of white spot syndrome virus genes [79]. Trypsin in the myocardium was able to trigger acute myocarditis following influenza A virus infection [80]. Overexpression and RNA silencing studies revealed that peroxidase was involved in the production of hepatitis C virus particles [81]. Moreover, RT-qPCR results indicated that the expression levels of all genes were down-regulated in BC9 and A35 following BmNPV infection (Fig 6). Therefore, we speculated that the down-regulation of these genes might be involved resistance to BmNPV infection.

Hypothesized modal analysis of the roles of the screened DEGs in silkworm resistance to BmNPV infection pathway
We hypothesized that the 22 DEGs discussed above played a role in the process of host response to BmNPV infection. The endocytosis process is triggered when the BmNPV nucleocapsid containing envelope binds to the cytomembrane. Vacuolar ATP synthase is activated by LPH to promote the fusion of the envelope and endosome thereby releasing the nucleocapsid into the cytoplasm. This process can be promoted by PAN1 and otoferlin. However, the transmembrane transport channel is an alternative pathway for virus to enter the cytoplasm, a process which can be facilitated by BAT1. The released nucleocapsid is transported into the nucleus with the help of the cytoskeleton (EFP). Once viral DNA is released into the nucleus, it will utilize host nucleotides to complete replication. In the final step of replication, viral DNA has to rely on host cell amino acids for assembly on the cytoskeleton (Table 5) [82]. In the cytoplasm, EFHP, ASGPA, ALT2, U2AF, ACT5 and ZRDP play an important role in facilitating virus replication, although the exact mechanism is still unclear. Moreover, transmembrane protein, MFS, is induced by BmNPV to increase cell volume, leading to lysis and cell death. We speculated that the down-regulation of these genes may affect the entry of virus into host cells and virus replication. The apoptosis process could also be triggered by PDK to inhibit BmNPV from further infecting other cells once BmNPV entered a host cell (Fig 7).
Taken together, our results provide useful information on silkworm resistance to BmNPV infection.
Supporting Information S1