extfieldsofstudy: Medicine Biology

s2fieldsofstudy: Biology Engineering

text:
Bacterial Abundance and Community Composition in Pond Water From Shrimp Aquaculture Systems With Different Stocking Densities

In shrimp aquaculture, farming systems are carefully managed to avoid rearing failure due to stress, disease, or mass mortality, and to achieve optimum shrimp production. However, little is known about how shrimp farming systems affect biogeochemical parameters and bacterial communities in rearing water, whether high stocking densities (intensive system) will increase the abundance of pathogenic bacteria. In this study, we characterized bacterial communities in shrimp ponds with different population densities. Water quality, such as physical parameters, inorganic nutrient concentrations, and cultivable heterotrophic bacterial abundances, including potential pathogenic Vibrio, were determined in moderate density/semi-intensive (40 post-larvae m-3) and high density/intensive shrimp ponds (90 post-larvae m-3), over the shrimp cultivation time. Free-living and particle-attached bacterial communities were characterized by amplicon sequencing of the 16S rRNA gene. Suspended particulate matter (SPM), salinity, chlorophyll a, pH, and dissolved oxygen differed significantly between semi-intensive and intensive systems. These variations contrasted with the equal abundance of cultivable heterotrophic bacteria and inorganic nutrient concentrations. Bacterial communities were dominated by Gammaproteobacteria, Alphaproteobacteria, Flavobacteriia, Bacilli, and Actinobacteria. Halomonas and Psychrobacter were the most dominant genera in the particle-attached fractions, while Salegentibacter, Sulfitobacter, and Halomonas were found in the free-living fractions of both systems. Redundancy analysis indicated that among the observed environmental parameters, salinity was best suited to explain patterns in the composition of both free-living and particle-attached bacterial communities (R2: 15.32 and 12.81%, respectively), although a large fraction remained unexplained. Based on 16S rRNA gene sequences, aggregated particles from intensive ponds loaded a higher proportion of Vibrio than particles from semi-intensive ponds. In individual ponds, sequence proportions of Vibrio and Halomonas displayed an inverse relationship that coincided with changes in pH. Our observations suggest that high pH-values may suppress Vibrio populations and eventually pathogenic Vibrio. Our study showed that high-density shrimp ponds had a higher prevalence of Vibrio, increased amounts of SPM, and higher phytoplankton abundances. To avoid rearing failure, these parameters have to be managed carefully, for example by providing adequate feed, maintaining pH level, and removing organic matter deposits regularly.

INTRODUCTION
Litopenaeus vannamei (Boone, 1931), white-leg shrimp, is an important species in aquaculture industry, which is widely reared in subtropical to tropical regions. They grow rapidly, have high survival rates even at high densities, possess a wide tolerance range of salinity and temperature, and can be cultivated in indoor (tanks or recirculating aquaculture systems) or outdoor facilities (ponds; Cuzon et al., 2004). Anderson et al. (2016) reported that world shrimp production in 2015 was 4.2 million tons of which L. vannamei contributed 75%. However, severe economic loss due to massive shrimp mortality has been occurring since 1987. For the period 1987-1994, economic loss reached USD 3 billion, representing 40% of the total production capacity of the industry (Brun et al., 2009;Walker and Winton, 2010). Moreover, a recent emerging bacterial disease, known as acute hepatopancreatic necrosis disease (AHPND), has occurred in Asian and Mexican shrimp aquaculture, causing an annual loss amounting to USD 1 billion with shrimp mortality exceeding 70% (Global Aquaculture Alliance [GAA], 2013, De Schryver et al., 2014Soto-Rodriguez et al., 2015).
The intensification of the shrimp industries resulted in changes of farming systems and sustainability. Until the year 2000, 70% of shrimp farming in Indonesia was conducted in extensive pond systems (Kautsky et al., 2000), but recently changed to semi-intensive or intensive systems (Kementerian Kelautan dan Perikanan [KKP], 2015). Different shrimp stocking densities affect rearing processes, including a defined nutritional input and shrimp production. Robertson et al. (1993) reported that the shrimp growth rate increased progressively as feeding frequency increased from 1 to 4 times per day. However, only a portion of the nutrients in the feed is consumed, assimilated, and retained as shrimp biomass. Shrimps only incorporate 24 to 37% of nitrogen and 11 to 20% of phosphorus from the feed into their bodies (Funge-Smith and Briggs, 1998;Sahu et al., 2013). In addition, 15% of nitrogen losses occurs during the first 2 h of immersion of feed pellets into the pond water (Smith et al., 2002). These unused nutrients will lead to a change of pH and dissolved oxygen (DO) in the water column and pond sediment, eutrophication, proliferation of bacteria and plankton, and an increase of particulate organic matter (Avnimelech et al., 1994;Martin et al., 1998). Furthermore, Kautsky et al. (2000) reported that the risk of shrimp diseases often increased with culture intensity and high stocking densities. Thus, shrimp diseases, along with bad water quality management, might threaten shrimp farming sustainability.
Free-living (FL) and particle-associated (PA) bacterial communities from the same water sample can be distinct (Bidle and Fletcher, 1995;Riemann and Winding, 2001;Jain and Krishnan, 2017;Yang et al., 2017). As organic-rich particles, aggregates provide a suitable habitat for microorganisms to take up nutrients, and shelter from predators, as well as from destructive physical factors (Lyons et al., 2010;Kramer et al., 2013). In addition, aggregates can accommodate higher bacterial abundance and diversity than the adjacent water column (Kiørboe et al., 2003). Due to water movement in ponds, organic matter agglomerates and forms large flocs or aggregates (Hargreaves, 2013;Avnimelech, 2014), which might facilitate bacterial settlement and proliferation. Previous studies about bacterial community compositions (BCCs) in shrimp pond waters have yet to analyze FL and PA fractions separately (Sombatjinda et al., 2011;Zhang D. et al., 2014;Xiong et al., 2016;Hou et al., 2017), resulting in a paucity of information on BCC in both fractions.
Although the causative agents for bacterial diseases in shrimps have been identified (Lavilla-Pitogo et al., 1998;Kharisma and Manan, 2012;Soto-Rodriguez et al., 2015;Xiao et al., 2017), preventive efforts to minimize disease outbreaks seem to still be ineffective. Studies in temperate ecosystems showed that pathogenic bacteria have been found in aggregates (Lyons et al., 2005;Froelich et al., 2013). Based on that evidence, we propose that particle abundance can be used to estimate the potential proliferation of pathogenic bacteria in shrimp farming, and that controlling aggregates may become an effective tool to manage the spread and survival of pathogens. Therefore, it is necessary to investigate water quality parameters, bacterial abundance, and BCC from different shrimp farming systems for both FL and PA fractions, over shrimp cultivation time.
The current study aims to comprehend the effect of shrimp farming systems of different intensity on water quality parameters and BCC, to elucidate factors affecting bacterial communities, and to evaluate the presence of pathogenic bacteria, including pathogenic V. parahaemolyticus in the FL and the PA/aggregates fraction. We hypothesize that (i) different shrimp farming systems affect pond water quality, including suspended particulate matter (SPM) loading, bacterial abundance and community composition, and (ii) SPM loads more bacterial cells, including pathogenic bacteria.

Sample Collection and Sampling Sites
Water samples were collected between 9-10 a.m. from no water-exchange and plastic lining ponds (square shape in size 2700-3000 m 2 , water depth 1.3-1.5 m) of semi-intensive (40 post-larvae m −3 , three ponds) and intensive (90 post-larvae m −3 , three ponds) systems, during a cycle of shrimp rearing at day 10, 20, 30, 40, 50, 60, and 70 (September to November 2016). Shrimp ponds were located in Rembang Regency, Central Java, Indonesia (−6 • 37 41.13 S 111 • 30 1 E and −6 • 42 11.66 S 111 • 21 54 E, for semi-intensive and intensive system, respectively). Two liters of water were taken from 5 points of each pond at 1 m depth, and then mixed. Two liters of the mixed water were prepared for bacterial abundance and community analysis; 3 liters for SPM and inorganic nutrient analysis. Samples were stored in a cold and dark container, and transported to the laboratory at Diponegoro University, Semarang, Indonesia, for immediate analysis. Remaining water was used for physical parameter measurement. Total harvest of each pond was recorded at the end of shrimp rearing.

Environmental Parameters
Salinity, temperature, pH, chlorophyll a, DO, and turbidity were measured ex situ, using calibrated Manta Eureka 2 multi-probes (Eureka Environmental Engineering, TX, United States).
Of each water sample, 0.3 L were filtered through a 0.45 µm syringe filter and the filtrate was poisoned with 1.2 mL of a 0.35 g L −1 HgCl 2 solution for inorganic nutrient analysis (ammonium: NH 4 + , nitrate: NO 3 − , nitrite: NO 2 , phosphate: PO 4 3− and silicate: SiO 4 4− ). The samples were stored at −20 • C until analysis. Inorganic nutrient measurements were done in triplicates at the Laboratory of Chemistry, Research Center for Oceanography (LIPI), Jakarta, Indonesia, according to the colorimetry method by Strickland and Parsons (1972), with a Shimadzu UV-1800 spectrophotometer.
Suspended particulate matter was measured as dry mass on pre-combusted GF/F filters (porosity 0.7 µm, ø 47 mm, VWR, France) in triplicates after filtration of a known volume of water sample (0.1-0.5 L). Weight of the filters was determined using a precision balance (ME 36S, Sartorius, Göttingen, Germany) after drying the filters for 24 h at 40 • C.
Samples for total bacterial cell counts were prepared by fixing 50 mL water with 4% v/v paraformaldehyde and stored at 4 • C for 24 h. The dilution factors for the samples were as follows: undiluted for day 10 samples, 5 × 10 −1 for day 20 and 30 samples, 10 −1 for day 40, 50, and 60 samples, and 5 × 10 −2 for day 70 samples. Ten milliliters of diluted fixed samples were subsequently filtered through 3.0 and 0.2 µm polycarbonate filters (ø 47 mm and 25 mm, respectively, Whatman, Dassel, Germany) to determine bacterial cells in the particle-attached (PA) and the free-living (FL) fractions, respectively. For the FL fractions at day 30, 40, and 70, only 5 ml filtrates were filtered. Filters were air dried and stored at −20 • C for further staining. A 4 ,6-diamidino-2-phenylindole (DAPI) staining was performed according to Kepner and Pratt (1994) for selected samples (10, 40, 50, 60, and 70 days samples). Filters were stained with 1 µg mL −1 of DAPI solution for 5 min, then washed in 80% ethanol and rinsed with sterile distillated water. Stained filters were air dried in the dark for 30 min, and then mounted with 10 µL of mounting solution consisting of 3:1 Citifluor AF mounting medium (Citifluor Ltd., London, United Kingdom) and Vectashield (Vector Laboratories Inc., Burlingame, United States). Bacterial cells, as well as size and number of aggregates were observed under a fluorescence microscope Axio Imager.D2 (Zeiss, Jena, Germany) at 1000× magnification. Bacterial cell abundance in FL filters was calculated from 30 photos per filter, using the free software ImageJ. Bacterial cells in PA filters were manually counted from 10 aggregates per filter of similar size. Sizes of aggregates were determined using a net micrometer grid (12.5 mm × 12.5 mm, divided into 10 × 10 fields, which is equal to 15,625 µm 2 at 1000× magnification). The cell average per filter was then divided by volume of filtered samples multiplied by dilution and factor of effective filter area (which was 31888 at magnification 1000×, for a filtration funnel with a diameter 25 mm, Millipore, Darmstadt, Germany) and number of aggregates (4 aggregate sizes per filter, Supplementary Table 1), for the FL and the PA fraction, respectively.

Molecular Analysis of Bacterial Communities
Five hundred milliliters of water samples were filtered subsequently through 3.0 and 0.2 µm polycarbonate filters (ø 47 mm, Whatman, Dassel, Germany) for PA and FL bacterial fractions, respectively. Genomic DNA was extracted according to Nercessian et al. (2005). DNA pellets were dissolved in TE buffer (10 mM Tris-HCl, 1 mM EDTA, pH 8.5). DNA concentrations were measured photometrically and checked for purity (ratio of light absorption at 260-280 nm) using a nanoquant plate reader (Infinite M200 Pro, Tecan, Germany). 16S rRNA gene amplification was performed from genomic DNA extracts from days 10, 40, 50, 60, and 70, considering cultivable bacterial abundance information (heterotrophic bacteria and potential pathogenic Vibrio) and bacterial disease evidence (white feces disease), which occurred at previous rearing cycles between 50 and 65 days of rearing (personal communication with shrimp farm owners).
DNA sequences of the V3-V4 hypervariable region of the 16S rRNA gene were obtained from amplicon sequencing with the primer set S-D-Bact-0314-b-S-17 (5 -CCTACGGGNGGCWGCAG-3 )/S-D-Bact-0785-a-A-21 (5 -G ACTACHVGGGTATCTAAKCC-3 ; Klindworth et al., 2013). Sequencing at LGC genomics (Berlin) was done on an Illumina MiSeq using the V3 Chemistry (Illumina) in a 2 × 300 bp paired-end run. Demultiplexing, i.e., grouping of sequences by sample, and the removal of the primer sequences from the raw paired-end reads were performed by LGC genomics (Berlin, Germany). Further bioinformatic analysis steps were carried out at ZMT, Bremen, Germany, according to Hassenrück et al. (2016). Sequences were quality trimmed with a sliding window of 4 bases and a minimum average quality of 15 with trimmomatic v.033 (Bolger et al., 2014). Quality trimmed sequences were merged using PEAR v0.9.8 . The swarming approach was used to cluster OTUs using swarm v2.1.1 (Mahé et al., 2014). For each OTU, one representative sequence (seed sequence) was taxonomically classified with SINA (SILVA Incremental Aligner) v1.2.11 using the SILVA rRNA project reference database (release 128) at a minimum alignment similarity and quality of 0.9 and a last common ancestor consensus of 0.7 (Pruesse et al., 2012). Unwanted lineages (such as archaea, chloroplasts, and mitochondria), as well as singletons and doubletons, i.e., OTUs occurring only once or twice in the data set, were removed. Rarefaction curves were calculated based on OTU richness and inverse Simpson index to assess the quality of sequence data sets. Samples with fewer than 500 sequences were excluded from the dataset for the analysis.
DNA sequence datasets were deposited in the European Nucleotide Archive (ENA) with the project accession number PRJEB26390, using the data brokerage service of the German Federation for Biological Data (GFBio, Diepenbroek et al., 2014), while biogeochemical data were archived in PANGAEA 1 .

Detection of Virulence Genes via PCR
Detection of toxin genes, i.e., transcriptional regulator (toxR), thermolabile haemolysin (tlh), thermostable direct haemolysin (tdh), Photorhabdus insect-related (pirA and pirB) was performed in a Mastercycler R (Eppendorf, Nexus gradient, Hamburg, Germany) using the set of primers described in Supplementary  Table 2. As these genes are common genes in several representative of the genus Vibrio (i.e., V. cholera, V. vulnificus, V. parahaemolyticus, V. alginolyticus, and V. owensii), we designed specific primer pairs which only amplify DNA sequences belonging to V. parahaemolyticus as this bacterium causes most of bacterial diseases in shrimp (Soto-Rodriguez et al., 2015;Heenatigala and Fernando, 2016;Xiao et al., 2017). Ten samples, which had a high amplicon abundance of Vibrio, were chosen for PCR analysis (Supplementary Table 3). The PCR conditions were as follows: a reaction mixture consisted of 2 µL (20 ng µL −1 ) of template, 2 µL of 10× PCR buffer B containing 15 mM MgCl 2 , 0.5 µL of 25 mM MgCl 2 , 0.5 µL of 0.2 mM dNTPs, 1 µL of 0.1 mM forward and reverse primer, 0.1 µL of 5U µL −1 Taq polymerase (all reagents provided by Roboklon EURx, Berlin, Germany) and 12.9 µL sterile distillated water. The amplification conditions were as follows: pre-denaturation at 95 • C for 3 min, followed by 40 amplification cycles consisting of denaturation at 95 • C for 30 s, annealing at 60 • C for 15 s, and extension at 72 • C for 30 s, and a final elongation at 72 • C for 5 min. Vibrio parahaemolyticus DSM 11058 (DSMZ, Braunschweig, Germany) was used as positive control for those toxin genes, while Vibrio vulnificus DSM 10143 (DSMZ, Braunschweig, Germany) served as negative control.

Data Analysis
To examine differences in environmental parameters and bacterial abundances between culturing intensity (intensive and semi-intensive systems), and among days, as well as the 1 https://doi.pangaea.de/10.1594/PANGAEA.889316 interaction between culturing system and sampling day, general linear mixed models (GLMM) were performed with shrimp pond as random factor. Data were log-transformed to achieve normal distribution prior to statistical testing. Tukey's post hoc tests were applied in cases where there were significant differences among sampling days and/or an interaction between sampling day and shrimp farming system. Forward model selection based on the Akaike information criterion (AIC), considering collinearity of the variables and variance inflation values, was used to determine environmental parameters, which best explained observed counts of cultivable heterotrophic bacteria and potential pathogenic Vibrio. Total shrimp harvest was tested with one-way ANOVA and total bacterial cell numbers from DAPI counting of FL and PA bacterial communities were analyzed using MANOVA. As first-level post hoc test, the effect of significant predictor variable (here: sampling day) in MANOVA of FL and PA bacterial numbers was then tested with individual GLMMs per fraction. Furthermore, pairwise comparisons between sampling days were conducted as second-level post hoc test. Repeated measures correlation between FL and PA bacterial cell numbers was estimated using rmcorr R package (Bakdash and Marusich, 2018).
Principal component analysis (PCA) was conducted to examine the relationship among environmental parameters, as well as cultivable bacterial abundances, and to characterize shrimp ponds of the different farming systems over time. Bray-Curtis dissimilarity coefficients were calculated to investigate variation of OTU compositions in FL and PA communities among ponds, as well as within and between the systems. Patterns of BCC in both fractions was visualized by non-metric multidimensional scaling (NMDS) based on Bray-Curtis dissimilarities using the function metaMDS of the vegan R package (Oksanen, 2017). Analysis of similarity (ANOSIM) was performed to test for differences in community composition in free-living and particle-attached fraction between systems over the time.
Redundancy analysis (RDA) was performed to explore relationships between bacterial communities and environmental variables in FL and PA fractions, for OTUs occurring in at least 10 samples, and with a sequence proportion of at least 1% in one sample. Data sets were separated by size fractions, and then centered log-ratio-transformed using aldex.clr (Fernandes et al., 2014) before being tested. Forward model selection based on the AIC, followed by collinearity and variance inflation inspection, was used to select the environmental parameters best suited to explain patterns in BCC. The significance of the individual parameters of the RDA models was assessed using restricted permutation tests to account for the repeated measurements within ponds. For the PA fraction, data from day 60 was excluded due to a missing value to ensure equal numbers of observation from each pond.
Population turnover, which occurred in the FL and PA fractions, was estimated for certain interesting OTUs of the genera Salegentibacter, Psychrobacter, Halomonas, and Vibrio, using the fraction of read abundance times cells method (FRAxC) proposed by Kevorkian et al. (2018), where bacterial cell numbers in the FL and PA of each selected sampling day were multiplied with the relative sequence proportion of 16S rRNA genes. We assumed that the biases associated with sampling procedure, DNA extraction, amplification, and cell counts did not vary over time. FRAxC results were also used to estimate correlation between estimated Vibrio abundance and selected environmental parameters (SPM, salinity, pH, temperature, ammonium, phosphate, and nitrite).

Shrimp Culturing Conditions
Total shrimp harvest differed significantly between semi-intensive and intensive systems (one-way ANOVA, F 1,4 = 111.43, p < 0.01). The intensive system showed a twofold higher production than semi-intensive system, which were 3,950 ± 284 kg and 1,990 ± 151 kg, for intensive and semi-intensive final harvest (mean ± SD, n = 3 per system). There was a high variability among the ponds within each system for turbidity, chlorophyll a, pH, salinity, DO, and inorganic nutrients during rearing time, which increased over time ( Table 1). Physical parameters such as salinity, pH, DO, chlorophyll a, and SPM were significantly different between the two systems (GLMM, p < 0.05, Supplementary Table 4). Among sampling days, we detected differences in temperature, salinity, pH, turbidity, chlorophyll a, SPM, cultivable heterotrophic bacteria and potential pathogenic Vibrio (GLMM, p < 0.05, Supplementary Table 4). The interaction between day and farming system was also significant for temperature, salinity, and turbidity. However, there were no significant differences between semi-intensive and intensive systems for dissolved inorganic nutrients, temperature, and turbidity (GLMM, p < 0.05, Supplementary Table 4).
In both systems, pH decreased rapidly after day 40, but when limestones were added into the ponds, it increased from 7.92 ± 0.13 to 8.38 ± 0.53 and 7.63 ± 0.46 to 8.11 ± 0.13, in semi-intensive and in intensive system, respectively. Salinity decreased from 35.37 ± 0.15 to 32.44 ± 3.26 PSU and 38.43 ± 0.78 to 33.57 ± 0.45 PSU, in semi-intensive and in intensive system, respectively, due to freshwater input. Overall, the intensive system showed significantly higher salinities, SPM, and chlorophyll a concentration than the semi-intensive system, as well as higher variability of inorganic nutrients ( Table 1).
Principal component analysis for environmental parameters and cultivable bacterial abundances was able to retrieve 42.61% of the variation among ponds on the first two principal components (Figure 1). According to water parameters and cultivable bacterial abundances, there was a clear separation by shrimp farming systems based on rearing time. At the beginning of rearing (days 10 and 20), ponds of intensive and semi-intensive systems were clustered with high values of the first principal component, which was driven mainly by high salinity. Afterward, they were separated by the second principal component with ponds of the intensive system characterized by high concentrations of chlorophyll a. The abundances of THB and TPPV highly correlated to SPM (Spearman correlation, rho: 0.87 and rho: 0.40 for THB and TPPV, respectively).
Aggregates containing bacterial cells were found in different sizes. At day 10, small aggregates appeared which were formed by bacterial cells and their exudates. After 40 days, aggregates were composed of plankton, bacterial cells and bacterial exudates, causing the increase of aggregate sizes in both systems, as well as bacterial cell numbers. Aggregates of sizes between 937.5 and 1,406.5 µm 2 were most abundant in both system over time, with the highest aggregate number at day 10 (511 ± 63 aggregates mL −1 ) and at day 40 (494 ± 27 aggregates mL −1 ), containing 86 ± 3 and 72 ± 6 bacterial cells per individual aggregate for semi-intensive and intensive pond waters, respectively (mean ± SD; n = 3; Supplementary Figure 1 and Supplementary  Table 1). Free-living (FL) bacterial cells of both systems increased and peaked at day 60 in concentration 3.8 × 10 7 cells mL −1 and 5.0 × 10 7 cells mL −1 , for semi-intensive and intensive systems, respectively. Total particle-attached (PA) bacterial cell numbers were steady after day 40, with 6.8 × 10 5 cells mL −1 and 6.4 × 10 5 cells mL −1 , for semi-intensive and intensive systems, respectively ( Table 2). The FL cell numbers were positively correlated to the PA cell numbers (repeated measure correlation, Spearman correlation, rho = 0.50, df = 23, p = 0.01) and differed only among days (Supplementary Tables 5, 6).

Bacterial Community Analysis
After removal of low-quality reads, 8,994-78,615 sequences, with on average 34,471 sequences per sample, were obtained from 59 samples. One sample was removed due to low sequence counts (<500 sequences). A total of 77,433 OTUs was obtained, ranging from 426 to 2465 OTUs, with an average of 1,312 OTUs per sample (Supplementary Figure 2).
Most abundant bacterial OTUs in both systems belonged to the classes Acidimicrobiia, Actinobacteria, Alphaproteobacteria, Bacilli, Cyanobacteria, Flavobacteriia, and Gammaproteobacteria (Figure 2). Among them, the genera Alteromonas, Erythrobacteraceae, Exiguobacterium, Halomonas, Vibrio, Pseudoalteromonas, Psychrobacter, Salegentibacter, and Sulfitobacter were present in every sample. Cyanobacteria, such as Synechococcus and Cobetia, were frequently present in the semi-intensive system. The highest sequence proportion in the FL fraction of both systems belonged to Salegentibacter,    Sulfitobacter, and Halomonas, while in the PA fraction Psychrobacter, Vibrio, and Halomonas were dominant, where the latter comprised up to 80% of the sequences in samples from the semi-intensive system after 40 and 50 days (Figure 2). In the SILVA reference database (SILVA 128 version), Halomonas OTU1 has a sequence similarity of 98 and 100% to H. aquamarina and H. meridiana, respectively. Free-living and particle-associated bacterial communities from the same water sample were very different from each other as indicated by Bray-Curtis dissimilarity coefficients (Table 3). In addition, there was high variability in BCC within the same fraction among replicate ponds, which lead to a wide range in Bray-Curtis dissimilarity coefficients (0.27-0.97 and 0.33-0.83 for PA and FL fractions, respectively) between intensive and semi-intensive systems ( Table 4).
NMDS showed a highly heterogeneous composition of the bacterial communities in the FL and the PA fraction for both systems (Figure 3). Analysis of similarity (ANOSIM) confirmed that there was no detectable pattern in FL or PA the BCC between system and day (Supplementary Table 7).
Redundancy analysis revealed that environmental variables explained 20.53 and 36.77% of the BCC in the FL and the PA fractions, respectively (Supplementary Figure 3). Among the observed environmental parameters, salinity was best suited to explain patterns in the composition of both FL and PA bacterial communities (R 2 > 10%). Furthermore, chlorophyll a and nitrate had a minor effect on BCC in the FL fractions ( Table 5).

Toxin Gene Assay, Vibrio Occurrence and Correlation to Environmental Parameters
All tested samples were negative for toxR, tlh, and tdh. These results indicate that no common pathogenic V. parahaemolyticus occurred in both systems, even though the given samples displayed high sequence proportions of Vibrio, specifically OTU4. The sequence of Vibrio OTU4 was identical to 16S gene sequences in the SILVA reference database belonging to V. mytili, V. diabolicus, and V. parahaemolyticus. The estimated abundance (FRAxC) of the dominant Vibrio OTU differed strikingly from that of the other dominant OTUs, such as Halomonas, Psychrobacter, and Salegentibacter, suggesting an inverse relationship. When FRAxCs of Vibrio were overly high, the proportions of the other genera were low. In the PA fraction, FRAxCs of Vibrio were negatively correlated to FRAxCs Halomonas (Figure 4).
FRAxCs of Vibrio in both fractions were positively correlated to SPM, temperature, ammonium, NP ratio, and TPPV, but they were negatively correlated to pH and salinity. In addition, Vibrio was negatively correlated to inorganic nitrite and phosphate in the FL fraction, although in the PA fraction the opposite trend was observed. In general, all correlation coefficients did not suggest strong correlation with values between −0.36 and 0.38 (Supplementary Table 8).

DISCUSSION
No water discharge shrimp aquaculture resulted in an accumulation of organic matter, through feed pellets, shrimp feces and other organic materials. Previous studies reported that degradation of organic matter, leaching of feed pellets, and release of ammonium from feces increase inorganic nutrients in the water column (Funge- Smith and Briggs, 1998;Smith et al., 2002;Burford and Lorenzen, 2004). These nutrients support bacterial and phytoplankton growth, which then leads to simultaneous change of physical parameters such as pH, DO, SPM, and turbidity (Burford and Williams, 2001;Burford et al., 2003).
FIGURE 2 | Contribution of the most abundant bacterial operational taxonomic units (OTUs) in the semi-intensive (light gray to black) and the intensive systems (yellow to dark red). Taxonomic affiliation for OTUs is provided for genus (Right) and class (Left) levels. S1-3 and T1-3: replicate pond for the semi-intensive and intensive ponds, respectively. Days of rearing are indicated below the pond symbols. FL, free-living fraction; PA, particle-attached fraction.
In this study, shrimp pond water parameters and respective BCC in intensive and semi-intensive systems were compared. Contrary to our expectations, higher shrimp density in intensive farming systems did not significantly increase inorganic nutrients in the water column. We suggest that excessive nutrients were rapidly taken up by bacteria and phytoplankton as their abundance increased more strongly in the intensive system. This was also reflected in high SPM values in the intensive system due to the aggregation of phytoplankton cells. Phytoplankton can take up total ammonia nitrogen (TAN), while heterotrophic bacteria perform nitrification and nitrogen assimilation (Avnimelech, 2014). In pond water, photosynthesis by phytoplankton and cyanobacteria increases DO, while respiration and microbial activities, such as nitrification and sulfur oxidation decrease pH (Boyd and Tucker, 2014). The higher number of phytoplankton (which was indicated by Chl a) in intensive ponds might increase DO concentration. We supposed that photosynthesis might increase oxygen concentration. However, we are aware that our DO measurements only constitute daylight conditions and that due to the higher organic matter loading of the intensive system and in the absence of light, oxygen depletion via respiration processes may lower DO concentration below the levels of S, semi-intensive system; T, intensive; NA * indicated a missing sample for particleattached fraction. the semi-intensive system. Even though the intensive system resulted in higher abundances of phytoplankton and cultivable heterotrophic bacteria, we could show that the concentrations of harmful inorganic nutrients, such as ammonium and nitrite, DO, and pH were still far from lethal values for shrimps, as mentioned in other studies (Schuler, 2008;Avnimelech, 2014;Furtado et al., 2015). We observed fluctuation of particle (aggregate or flocs) number over rearing time. We supposed that lower pellet input, sinking or consumption of flocs might cause the fluctuation of particle numbers. Water movement and exudate secreted either by bacteria or plankton lead to particle agglomeration (Grossart et al., 2006;Gärdes et al., 2011), generating bigger aggregates/flocs. When the aggregates sink to the bottom, particle numbers might decrease, resulting low turbidity/clearer water. Unfortunately, we did not measure daily feed input, sinking rate of aggregates and particle consumption by shrimp.
Environmental factors shape the structure and function of microbial communities (Allison and Martiny, 2008). Previous studies have shown that the succession of microbial communities was influenced by combinations of chlorophyll a, total nitrogen (TN), PO 4 3− , C/N ratio , total phosphate, chemical oxygen demand (Zhang D. et al., 2014), and feed sources added into ponds (Qin et al., 2016). Of the parameters measured in this study, salinity was the most determinant variable which shaped BCC in the FL and particle-attached PA fractions, which concurred with previous studies (Yang et al., 2016;Hou et al., 2017;Kirchmann et al., 2017). We observed salinity fluctuation in both systems and supposed that rainfall and addition of sea water from reservoir pond to maintain rearing pond water level (∼140-150 cm) decreased salinity. During our study, rainfall could increase 10-20 cm of water level, while sea water was added more frequent in intensive ponds due to regular mud removal. However, it should be considered that the influences of salinity on the BCC in the FL fraction do not completely exclude the influence of other water parameters, such as unmeasured sulfur compound in water column like hydrogen sulfide (H 2 S) or sulfite which may explain the continually presence of Sulfitobacter in FL fraction.
Bacterial communities in intensive and semi-intensive systems were dominated by heterotrophic halophilic bacterial genera Alteromonas, Erythrobacter, Halomonas, Pseudoalteromonas, Salegentibacter, Sulfitobacter, and Psychrobacter. Halomonas and Psychrobacter, the salt-tolerant heterotrophic nitrifying bacteria (Chankaew et al., 2017), can oxidize ammonia in high concentration under saline conditions (Sangnoi et al., 2016). Sulfitobacter can oxidize sulfur and degrade high molecular weight dissolved organic matter (Bourne et al., 2004;Sosa et al., 2017). Other bacterial genera, such as Salegentibacter, Exiguobacterium, and Erythrobacter, were associated with shrimp shell degradation after molting (Sombatjinda et al., 2011). In our study, BCC differed from those reported by Sombatjinda et al. (2011) who analyzed L. vannamei earthen pond water with salinities close to freshwater. They identified Nitrosomonas, Flavobacterium, Exiguobacterium, Synechococcus, Burkholderia, Nitrosospira, and Nitrobacter as dominant bacteria (Sombatjinda et al., 2011). We assume that salinity caused the differences in BCC observed in our study, considering the salinity range of 32-38 PSU compared to 2-10 PSU in previous study (Sombatjinda et al., 2011). This is supported by a study of Chankaew et al. (2016), where common nitrifying bacteria, including Nitrosomonas, Nitrosococcus, Nitrosolobus, Nitrospira, Nitrococcus, and Nitrobacter, were absent in shrimp pond water with 20 PSU. Furthermore, the regular addition of commercial probiotics (personal communication with pond owners), which contained Bacillus sp., Pseudomonas sp, Nitrosomonas sp., Aerobacter sp., Nitrobacter sp., and Nitrosococcus sp. did not affect BCC in the FL and PA fractions of both systems, indicating FIGURE 3 | Non-metric multidimensional scaling (NMDS) plot of bacterial community composition (BCC) in the free-living (FL) and particle-attached (PA) fraction. Information on environmental parameters was added to the NMDS plot using envfit. that nitrification processes were generated by other microbes. We propose that in our study the dominance of particular heterotrophic halophilic bacteria may result in higher uptake of inorganic nutrient such as ammonium and nitrite (Sangnoi et al., 2016). Therefore, even in the intensive pond waters inorganic nutrients remained low. In our study, BCC changed over time, following a "resilience scenario" (Allison and Martiny, 2008), in which a replacement of bacterial taxa occurred due to environmental change, followed by a quick return to its pre-disturbance composition. This phenomenon can be clearly seen in the particle-attached fraction of the intensive system. When pH decreased below 8 at day 40 and 50, Vibrio replaced Halomonas as the most abundant OTU, but the sequence proportions of Halomonas recovered, after the pH increased at day 60 and 70. In contrast, BCC in the FL and the PA fractions of the semi-intensive system during the same period was similar while experiencing stable pH (above 8). We propose that pH was a disturbing factor for heterotrophic halophilic bacteria in our shrimp systems and lead to the change of BCC only in the intensive system. This is supported by Krause et al. (2012) who reported that slight changes in pH caused compositional shifts in marine bacterial communities. Vibrio, a potential opportunistic pathogen for L. vannamei (Heenatigala and Fernando, 2016;Rungrassamee et al., 2016), has been found in higher proportions in the PA fractions of the intensive system. This might indicate that the particulate fraction, specifically marine aggregates, can accumulate potentially pathogenic bacteria, as suggested by other studies (Lyons et al., 2005;Froelich et al., 2013). Our hypothesis that the intensive TABLE 5 | Contribution and significance of observe environmental factors which explain the variation in community composition of the free-living (FL) and the particle-attached (PA) fractions based on redundancy analysis (RDA).

Source of variation
Adjusted R  system seemed to be more vulnerable to Vibrio outbreaks is supported by high abundance and recurrent presence of Vibrio spp. only in the intensive system. Therefore, it is necessary to maintain SPM and aggregate abundance and to avoid massive Vibrio growth due to the fact that Vibrio can convert from non-virulent to virulent under certain cell density threshold or if dramatic environmental changes occur (Zhou et al., 2012). Nevertheless, toxin genes (toxR, tlh, and tdh) of V. parahaemolyticus (Makino et al., 2003), as well as pirAB toxins for AHPND (Xiao et al., 2017) were not present in this study. We argue that other halophilic, potentially probiotic, bacteria might inhibit chromosome II replication, where those genes are located (Makino et al., 2003), as proposed by Defoirdt et al. (2011). As no other known probiotic taxa were detected in this study, we hypothesize that the presence of Halomonas in both systems might inhibit Vibrio. Recently, Halomonas aquamarina has been applied as probiotic in L. vannamei culture to oxidize ammonium and to prevent Vibrio growth, which then leads to an increase of survival rates of L. vannamei (Zhang et al., 2009;Suantika et al., 2013;Sangnoi et al., 2016). However, other unobserved microorganisms, such as bacteriophage, Bdellovibrio, Saccharomyces, Streptococcus, Streptomyces, and protists might also suppress Vibrio (Munro et al., 2003;Alagappan et al., 2010;Chavez-Dozal et al., 2013;Lakshmi et al., 2013;Kongrueng et al., 2017). Considering final shrimp harvests, the intensive system is the more economically promising system. However, shrimp pond management needs be improved to maintain water quality as well as beneficial bacterial communities. Regular feed control and mud discard were likely to have affected BCC in our study. Without regular mud discard, sludge and organic matter degradation might increase oxygen demand, and as a result oxygen might become depleted in the bottom of the shrimp pond. Under these conditions, anaerobic degradation of organic matter produces H 2 S, and increases other toxic inorganic compounds such as ammonia and nitrite. Furthermore, sludge deposit in bottom pond might enlarge the anaerobic area which may reduce the habitable space for the shrimp. While strictly adjusted feeding in semi-intensive system seemed to lead to a more stable BCC, regular mud discard in the intensive system, which was conducted once to twice per day, might have prevented the accumulation of toxic compounds (i.e., hydrogen sulfide, ammonium, and nitrite), as indicated by the lower proportion of Sulfitobacter. We suggest that giving adequate feed, discarding sludge, maintaining pH, and salinity, might avoid abrupt water quality change, eliminate large anaerobic areas at the pond bottom, and may minimize proliferation of potential pathogens, especially Vibrio. Some aspects of shrimp pond management related studies are given in Hopkins et al. (1993); Funge- Smith and Briggs (1998), and Avnimelech and Ritvo (2003). Consequently, it is necessary for pond aquaculture to provide sufficient seawater and sludge reservoirs to avoid environmental degradation risks due to sludge removal. In view of BCC in FL and PA fraction over rearing period and the absence of virulence genes, we propose to apply "bioflocs" in shrimp aquaculture. This can be generated by adding carbohydrates such as molasses, rice bran or tapioca. Bioflocs might become an alternative food source for shrimps and be able to maintain ammonia and nitrite concentration in shrimp pond (Avnimelech, 2014). Moreover, because bacteria from commercial probiotics were not detected in any samples, we suggest that the regular application of these probiotics in such high salinity rearing pond is not necessary to be done. At the end, operational shrimp rearing costs, especially for feed pellet and commercial probiotics, can be reduced.

CONCLUSION
Different stocking densities influenced water quality parameters, especially SPM, DO, chlorophyll a, and pH. High salinity shaped the BCC in the PA and FL fraction, favoring heterotrophic halophilic bacteria which facilitate an optimum uptake of inorganic nutrients and prevent increase of Vibrio growth. Despite a wide variability of BCCs in intensive and semi-intensive systems, the abundance of Vibrio was higher in the intensive system, which may therefore be more vulnerable to disease outbreaks. Monitoring BCC, especially the PA fraction's, and larger particulates such as aggregates in shrimp pond waters may potentially prevent disease outbreaks such as Vibriosis, white feces disease, white-tale disease, and AHPND.

ETHIC STATEMENT
This study did not sacrifice any organism, nor introduce viable pathogenic bacteria to the environment. All sampling procedures were already informed to and have been agreed upon by shrimp pond owners.

AUTHOR CONTRIBUTIONS
YA, JH, and AG designed the project. YA collected the samples and conducted the in situ measurement with the logistic support of AT and AK. YA and CH completed the statistical analysis. YA prepared the manuscript with input from all co-authors.

FUNDING
We acknowledge financial support from the Leibniz Association. YA gets a funding for pursuing his Ph.D. with the RISET-Pro scholarship from Indonesian Ministry of Research, Technology and Higher Education.


extfieldsofstudy: Biology Medicine

s2fieldsofstudy: Biology Medicine

text:
Inhibition of long non-coding RNA HOTAIR enhances radiosensitivity via regulating autophagy in pancreatic cancer

Background Resistance to radiation therapy is still a challenge for treatment of pancreatic cancer(PC). Long non-coding RNAs (lncRNA) HOTAIR has been found to play a oncogenic role in several cancers. However, the correlation between HOTAIR and radiotherapy in PC is still unclear. Methods TCGA data was collected to analyze the expression of HOTAIR and its relationship with PC progression. A series of functional experiments were conducted to explore the role of HOTAIR in PC radiosensitivity and its underlying molecular mechanisms. Results By the analysis of the TCGA data, we found HOTAIR expression in PC tissues was significantly higher than normal tissues and associated with tumor progression. The function analysis showed HOTAIR was enriched in biological regulation and response to stimulus. And in vitro study, the expression of HOTAIR was increased in PANC-1 and AsPC-1 cells after radiation. We identified that HOTAIR knockdown could enhance radiosensitivity and influence autophagy by up-regulating ATG7 expression in PC cells. By futher rescue experiments using rapamycin, activation of autophagy could reversed the the inhibition of cell proliferation and colony formation, as well as promotion of apoptosis mediated by HOTAIR knockdown, indicating that HOTAIR knockdown promoted radiosensitivity of PC cells by regulating autophagy. Conclusion Our finding revealed the the regulatory role of HOTAIR in radiosensitivity and provided a a new sight to improve radiotherapy effciency in PC.

Introduction
Pancreatic cancer (PC) is a kind of highly invasive potential tumor. The survival rate of PC patients who have not performed radical resection is <2%, which is only 24% even after early radical resection. 1 PC patients who have not undergone radical resection are with a median survival time <6 months. In the western countries, the death caused by PC is the fourth, after lung cancer, colorectal cancer, and breast cancer. 2 Radioresistance is one of the main reasons for poor prognosis of PC patients. How to improve the sensitivity of PC with radiation is a clinical problem to be solved.
Long non-coding RNA (lncRNA) is of >200 nt in length, which is located in the nucleus and cytoplasm. 3 LncRNA has stability splice site, and it is with a specific secondary structure and expression of the space-time specialty. It can regulate gene expression in epigenetics, transcription, and post-transcription level. 3 Moreover, some studies indicated that lncRNA is related to cancer radioresistance, for example, lncRNA NKLA can increase laryngeal cancer cell radioresistance. 4 lncRNA HOTAIR is related to gynecologic cancer cell radioresistance. 5 lncRNA UCA1 can increase prostate cancer radioresistance. 6 A num-ber of studies 7,8 about PC have demonstrated that lncRNA can promote PC genesis and development, such as lncRNA UCA1, lncRNA DUXAP10, and lncRNA HOTAIR. However, the functions of lncRNA in PC radioresistance are still unrevealed.
In addition, autophagy plays a vital role in cancer radioresistance. Cancer radiation, as a physical agent, can kill cancer cells. When tumor cells obtained high-energy radiation, their genetic material (DNA) is destroyed, and the ability of selfrenewal, division, and proliferation is reduced. 9 Cancer cells adapted these negative factors through autophagy to remain stable. 10 Autophagy is a cellular metabolism with nutrient deficiency or other metabolic stress. Maintaining degradation of intracellular metabolic balance or aggregation of proteins and organelles is the main function of autophagy. Autophagy can make cancer cells survive from radiation therapy. Studies have demonstrated that inhibition of autophagy in cancer cells, such as breast cancer cells, nasopharyngeal carcinoma cells, esophageal carcinoma cells, and lung cancer cells, increases the radiation sensitivity of tumor cells. [10][11][12][13] Studies have shown that in PC, autophagy can promote the development of the occurrence of PC. 14,15 However, the role of autophagy in PC is rarely reported. In the process of autophagy in cancer cells, lncRNA also plays an important role. Xiao et al found that lncRNA MALAT1 can increase the drug resistance of hepatic cellular cancer (HCC) cells by regulating autophagy through HIF-2α-MALAT1-miR-216b axis. 16 Another research found that lncRNA HILC may regulate apoptosis and invasion of gastric cancer cells by regulating autophagy. 17 To sum up, lncRNA may have an effect on PC cell radioresistance through regulating autophagy and PC cells. In the present study, we intended to find the relationship between autophagy and PC cell radioresistance, which may provide theoretical foundation for the clinical radiotherapy of PC.

Materials and methods
The Cancer genome atlas (TCga) analysis The pancreatic adenocarcinoma (PAAD) genomic data were downloaded from TCGA database (https://cancergenome.nih. gov). The expression data were pretreated by Bioconductor package. A total of 179 PC cases with detailed clinical information and the corresponding HOTAIR expression data were included in this study. The PC cases were divided into two groups according to the clinical information to analyze the relationship between HOTAIR expression and clinicopathological features. Kaplan-Meier method was used to draw the survival curves compared with the log-rank test.

Function analysis
The function enrichment analysis was conducted on the high and low HOTAIR expression in PC cases. The biological process category was used in the analysis. Statistical significance was considered as P<0.05.

Cell culture and transfection
The human PC cell lines PANC-1 and AsPC-1 cells were purchased from the Cell Bank of the Chinese Academy of Sciences (Shanghai, China). Cells were cultured in DMEM (Sigma-Aldrich Co., St Louis, MO, USA) with 10% FBS at 37ºc in humidified air with 5% CO 2 . All the vector plasmids such as pcDNA3.1-HOTAIR, si-HOTAIR, and negative control were chemically synthesized by Genechem Co. Ltd. (Shanghai, China). Cell transfection was performed using Lipofectamine 2000 (Thermo Fisher Scientific, Waltham, MA, USA).

Radiation treatment of cells
Cells were irradiated by Varian Clinac 23EX with 6-MV X-ray at 400 cGy/min dose rate. Cells were exposed to different doses of irradiation treatment (0, 2, 4, 6, and 8 Gy). After 24-96 hours, the cells were used for other analysis. For time-course experiment, the cells received 4 Gy and then were collected for quantitative real-time PCR (qRT-PCR) analysis every 2 hours for 24 hours after irradiation.

Proliferation assay
Cells were inoculated into a 96-well plate and transfected with si-HOTAIR or si-NC and then exposed to 4 Gy dose X-ray. At 0, 24, 48, 72, and 96 hours after radiation, 20 µL MTT reagent (Sigma-Aldrich Co.) was added into each well and incubated for another 4 hours at 37°C in the dark. The

Flow cytometric analysis of apoptosis
PANC-1 and AsPC-1 cells were inoculated into a six-well plate and transfected with si-HOTAIR or si-NC for 24 hours. Then the cells were treated with a 4 Gy dose X-ray. At 24 hours after irradiation treatment, annexin-V, fluorescein isothiocyanate, and propidium iodide apoptosis detection kit was used to measure cell apoptosis rate. Each experiment was conducted at least three times.

Colony formation assay
PANC-1 and AsPC-1 cells (1,000 cells/plate) transfected with si-HOTAIR or si-NC were seeded in six-well plates. At 24 hours after incubation, the cells were treated with different doses (0, 2, 4, 6, and 8 Gy) of X-ray irradiation at a dose rate of 3.5 Gy/min and seeded into 60 mm dishes until visible colonies appeared. The colonies were then stained with crystal violet (Santa Cruz Biotechnology Inc., Dallas TX, USA) and were counted.

autophagy-related detection
The autophagy flux detection was the same as that in our previous study. 18,19 statistical analysis Data analyses were performed using the SPSS 17.0 software (SPSS Inc., Chicago, IL, USA). The results are presented as mean ± SD; statistical analyses were performed using Student's t-test or analysis of variance. P<0.05 was considered statistically significant.

Expression profile and function analysis of hOTaiR in PC
To explore the expression profile of lncRNA HOTAIR in PC tissues, we analyzed the PAAD genomic data which can be downloaded from TCGA database (https://cancergenome.nih. gov). By analysis of the RNA-seq data, we found that HOTAIR was significantly upregulated in PC tissue samples than normal pancreatic tissues, where the raw data were obtained from GTEx ( Figure 1A). To analyze the relationship between

5264
Wu et al HOTAIR expression and clinicopathological features, the results showed that HOTAIR expression was correlated with tumor T stage and lymph node metastasis. The expression of HOTAIR was significantly higher in lymph node metastasis group ( Figure 1B) and was correlated with T stage ( Figure  1C). In addition, Kaplan-Meier analysis of 179 PC patients presented that PC patients with high HOTAIR expression has a poor prognosis than patients with relatively low HOTAIR expression ( Figure 1D). Function analysis of biological process for HOTAIR in PC showed that biological regulation and metabolic response to stimulus were the most enriched processes ( Figure 1E). Taken together, these results indicated that HOTAIR plays an important role in PC progression.
hOTaiR was upregulated after radiotherapy in PC samples and cells Growing evidence has shown that lncRNA can participate in the process of radiosensitivity. To explore the role of HOTAIR in radiotherapy, we performed RNA-seq analysis in PC samples which had received radiation therapy. There are 119 PC patients who had received radiotherapy and 44 PC patients without radiotherapy. Differential expression analysis revealed a set of differentially expressed genes, in which red points represented upregulated genes and green points represented downregulated genes ( Figure 2A). The analysis revealed HOTAIR was one of the significant upregulated gene expression ( Figure 2B). In in vitro experiments, we detected whether radiation treatment could influence the expression of HOTAIR in PC cells.
When treated with 4 Gy X-ray in PANC-1 and AsPC-1 cells, the expression of HOTAIR was significantly increased after 24 hours ( Figure 2C and D). Additionally, HOTAIR expression was measured in PANC-1 and AsPC-1 cells after different doses (0, 2, 4, 6 Gy) of radiation treatment for 24 hours ( Figure 2E and F). Taken together, these results indicated that the expression of HOTAIR in PC cells can be modulated by radiation treatment in time-and dose-dependent manners.

hOTaiR knockdown enhances radiosensitivity of PC cells
According to the upregulation of HOTAIR in response to radiation treatment, we next explored whether HOTAIR could affect the radiosensitivity in PC cells. First, we used three different siRNAs to knock down HOTAIR expression in PC cells. The si-HOTAIR3 has the highest transfection efficiency ( Figure 3A and B). Consequently, we used si-HOTAIR3 in the follow-up experiments. PANC-1 and AsPC-1 cells were exposed to 4 Gy radiation treatment transfected with si-NC or si-HOTAIR. The results indicated that knockdown of HOTAIR could further inhibit the proliferation of PC cells compared with only radiation treatment ( Figure 3C and D). Additionally, knockdown of HOTAIR or radiation therapy could induce apoptosis of PC cells. Combination of si-HOTAIR and radiation therapy significantly increased the apoptosis rates of PANC-1 and AsPC-1 cells compared to si-HOTAIR or 4 Gy+ si-HOTAIR treatment ( Figure 3E and F). Moreover, after exposed to various doses of X-ray, colony formation assay was performed in PC cells transfected with si-NC or si-HOTAIR. The results showed that the PC cell knockdown of HOTAIR were more sensitive to X-ray irradiation than si-NC group ( Figure 3G and H). These results indicated that HOTAIR knockdown can enhance radiosensitivity of PC cells.

hOTaiR can promote autophagy of PC cells
Growing evidence has demonstrated that autophagy is involved in radiosensitivity in cancer cells. And several studies have found that lncRNA can participate in the regulation of autophagy. We hypothesized that HOTAIR affected radiosensitivity through regulating autophagy in PC cells. To probe for the relationship between HOTAIR and autophagy, we detected the autophagy flux level when PC cells were overexpression of HOTAIR and knockdown of HOTAIR. The results showed that the ratio of LC3-II/LC3-I was increased. Because autophagy is a dynamic process, both induction of autophagosome formation and inhibition of autophagosomes degradation can lead to the accumulation of LC3-II. So we used bafilomycin to inhibit the degradation of LC3-II. The result showed that accumulation of LC3-II was enhanced in treated with bafilomycin and HOTAIR overexpression cells, which suggests that HOTAIR promoted the autophagosome formation ( Figure 4A and B). Then we utilized the LC3-fused GFP construct (pGFP-LC3) to investigate the distribution of autophagic marker protein LC3. We found that the GFP-LC3 puncta increased in the HOTAIR overexpression group ( Figure 4C and D). In addition, we observed the formation of autophagosome and autophagolysosome in PC cells transfected with pcDNA3.1-HOTAIR by the TEM method ( Figure 4E). HOTAIR markedly increased the number of autophagosomes ( Figure 4F). These results support that HOTAIR can promote autophagy flux of PC cells.
hOTaiR promotes autophagy by upregulating atg7 expression To explore the molecular mechanism by which HOTAIR promotes autophagy in PC cells, we used qRT-PCR to detect the main autophagy-related gene expression change in PC cells, which overexpressed HOTAIR. Among autophagyrelated genes of ATG3, ATG5, ATG7, ATG9, ATG10, ATG12, ATG14, and Beclin1, the expression of ATG7 was markedly increased by HOTAIR at mRNA level ( Figure 5A and B). In addition, the protein level of ATG7 was also upregulated in our immunoblot assay ( Figure 5C and D). On the contrary, knockdown of HOTAIR decreased ATG7 expression ( Figure  5E and F). ATG7 is a key regulator of the autophagosome formation and responsible for vesicle progression. 20 These

hOTaiR affects radiosensitivity of PC cells by regulating autophagy
Based on the results of HOTAIR regulating autophagy in PC cells, we further validated this regulatory axis between HOTAIR, autophagy, and radiosensitivity. The result showed that rapamycin treatment restored the autophagy levels in HOTAIR knockdown cells ( Figure 6A and B). PANC-1 and AsPC-1 cells were exposed to 4 Gy radiation transfected with si-NC or si-HOTAIR or combined si-HOTAIR and rapamycin (an activator of autophagy). MTT assay results showed that HOTAIR knockdown markedly decreased cell proliferation but rapamycin significantly reversed this effect in PC cells ( Figure 6C and D). In addition, rapamycin-mediated activation of autophagy in part reversed the apoptosis effect of si-HOTAIR ( Figure 6E and F). Knockdown of HOTAIR notably reduced colony formation rates in PANC-1 and AsPC-1 cells, which was reversed by rapamycin ( Figure 6G and H). These The effect of hOTaiR knockdown on radiosensitivity in PC cells. Notes: (A, B) The transfection efficiency of si-HOTAIR1, si-HOTAIR2, and si-HOTAIR3 was performed by qRT-PCR in PANC-1 and AsPC-1 cells. (C, D) MTT assay was performed to evaluate cell viability at 24, 48, 72, and 96 hours after PC cells transfected with si-hOTaiR or si-nC exposed to 4 gy radiation. (E, F) apoptosis rate was detected in PC cells transfected with si-hOTaiR or si-nC exposed to 4 gy radiation after 24 hours. (G, H) Colony formation assay was performed in PC cells transfected with si-nC or si-hOTaiR exposed to various doses of X-ray. *P<0.05.

Discussion
Accumulated evidence has identified that lncRNA can affect a number of biological processes of cancer cells, such as proliferation, migration, metastasis, apoptosis, and so on. 21,22 And more and more lncRNAs have been found to be aberrantly expressed in PC cells and involved in various steps of PC progression. Many lncRNAs have the potential value in diagnosis and prognosis in PC. 23 By next-generation sequenc-ing method, Müller et al 24 have found that 43 lncRNAs were abnormally expressed in PC tissues. In addition, Wang et al 25 aimed to probe for the potential biomarkers in 7,419 lncRNAs by microarray, in which chip platform was designed for screening lncRNA. The results showed that HOTTIP-005, XLOC_006390, and RP11-567G11.1 were significantly increased in PC tissue. HOTAIR was a lncRNA transcribed from the antisense strand of the HOXC gene locus and has been demonstrated to be overexpressed in various carcinomas. 26 In the present study, we identified HOTAIR expression profile in PC RNA-seq data from TCGA database. The advantage of TCGA database is its large sample size and detailed clinical parameter information. Therefore, the result is more persuasive. Our analysis result revealed that HOTAIR was upregulated in PC tissues and significantly related to prognosis. Besides, the function analysis showed that the abnormal expression of HOTAIR was associated with biological regulation and metabolic response to stimulus. Radiotherapy is one of the main treatment approaches for PC patients. 27 However, radioresistance is still the main limitation for the radiation treatment. Several lncRNAs have been found to play an important role in radiosensitivity. [28][29][30] The expression of lncRNA PVT1 has been identified to be increased after radiation in non-small-cell lung cancer (NSCLC) tissues. And PVT1 knockdown enhanced radiosensitivity of NSCLC cells by sponging miR-195. Shen et al 30 demonstrated that LincRNA-p21 knockdown enhances radiosensitivity of hypoxic tumor cells through HIF-1/Akt/ mTOR/P70S6K pathway. Yang et al found that ectopic expression of lncRNA GACAT3 can enhance the effect of radiotherapy in A549 cells. 28 These evidence revealed that lncRNA can influence radiosensitivity of cancer cells. From the result of our function analysis, we suspected that HOTAIR may participate in the biological regulation and response to stimulus in PC cells. In in vitro study, our results showed that HOTAIR knockdown enhances radiosensitivity of PC cells. This finding indicated that HOTAIR maybe a potential target for improving the efficiency of radiotherapy.
Next, our aim is to research the mechanism of HOTAIR regulating radiosensitivity in PC cells. A growing number of studies have found that autophagy is related to radiosensitivity. 31  Autophagy is a conserved lysosome-mediated and dynamic process that can generate energy and nutrients for maintaining cellular homeostasis. 32 It is widely accepted that autophagy plays a protective role in cancer cells during stressors. 33 Not only that, many studies have identified that autophagy can protect cells from IR-induced cellular damage and enhance radioresistance. 34,35 Inhibition of autophagy is important for the improvement of efficacy of radiotherapy in cancer. Recently, some studies found that lncRNA could activate autophagy by regulating Atg gene expression in cancer cells. 19,36 Based on these findings, we suspected that whether HOATIR as a lncRNA could influence radiosensitivity by regulating autophagy in PC cells. Our study found that the ratio of LC3-II/LC3-I and the formation of autophagosome were increased when HOTAIR was overexpressed, indicating that HOTAIR can promote autophagy in PC cells. Furthermore, we utilized the activator of autophagy in the rescue experiments. Our results showed that rapamycin could reverse the inhibition of cell proliferation and colony formation, as well as promotion of apoptosis in vitro mediated by HOTAIR knockdown. It revealed that HOTAIR knockdown promoted radiosensitivity of PC cells by regulating autophagy. Our finding was consistent with a previous study about non-coding RNA and radiotherapy, for example, Peng et al found that inhibition of miR-23b increases autophagy to promote radioresistance. 37 However, the detail molecular mechanism about this phenomenon should be further explored.

Conclusion
In conclusion, we identified that HOTAIR plays a key role in PC progression and it is correlated with radiotherapy through the analysis of the TCGA database. From in vitro approaches, we found that HOTAIR expression can be upregulated by radiation. And our study showed that HOTAIR knockdown can enhance radiosensitivity via regulating autophagy in PC cells. Our investigation of the relationship between HOTAIR and radiotherapy provided a novel highlight for lncRNA in cancer therapy.


text:
Next Generation Sequencing-Based Molecular Marker Development: A Case Study in Betula Alnoides

Betula alnoides is a fast-growing valuable indigenous tree species with multiple uses in the tropical and warm subtropical regions in South-East Asia and southern China. It has been proved to be tetraploid in most parts of its distribution in China. In the present study, next generation sequencing (NGS) technology was applied to develop numerous SSR markers for B. alnoides, and 64,376 contig sequences of 106,452 clean reads containing 164,357 candidate SSR loci were obtained. Among the derived SSR repeats, mono-nucleotide was the main type (77.05%), followed by di- (10.18%), tetra- (6.12%), tri- (3.56%), penta- (2.14%) and hexa-nucleotide (0.95%). The short nucleotide sequence repeats accounted for 90.79%. Among the 291 repeat motifs, AG/CT (46.33%) and AT/AT (44.15%) were the most common di-nucleotide repeats, while AAT/ATT (48.98%) was the most common tri-nucleotide repeats. A total of 2549 primer sets were designed from the identified putative SSR regions of which 900 were randomly selected for evaluation of amplification successfulness and detection of polymorphism if amplified successfully. Three hundred and ten polymorphic markers were obtained through testing with 24 individuals from B. alnoides natural forest in Jingxi County, Guangxi, China. The number of alleles (NA) of each marker ranged from 2 to 19 with a mean of 5.14. The observed (HO) and expected (HE) heterozygosities varied from 0.04 to 1.00 and 0.04 to 0.92 with their means being 0.64 and 0.57, respectively. Shannon-Wiener diversity index (I) ranged from 0.10 to 2.68 with a mean of 1.12. Cross-species transferability was further examined for 96 pairs of SSR primers randomly selected, and it was found that 48.96–84.38% of the primer pairs could successfully amplify each of six related Betula species. The obtained SSR markers can be used to study population genetics and molecular marker assisted breeding, particularly genome-wide association study of these species in the future.

Introduction
Betula alnoides Buch. Ham. ex D. Don (Betulaceae) is a fast-growing valuable indigenous tree species with multiple uses in the tropical and warm subtropical regions in South-East Asia and southern China. In China, its main natural distribution is in Yunnan, Guangxi, and Guizhou provinces [1]. The wood has beautiful texture, moderate density, low crack and deformation rate, and easy processing, making it suitable for floor timber and high-grade furniture making, high-class interior

Distribution of SSR in the Genome
A total of 64,376 contig sequences, about 60.47%, of 106,452 clean reads, contained SSR loci. Among these, 25,431 sequences contained only one SSR locus, and the other 38,945 sequences contained at least two SSR loci. In total 164,357 candidate SSR loci were obtained ( Table 1). The total size of examined sequences was 106,901,537 bp, and the size of the SSR repeat lengths ranged from 10 bp to 75 bp with a mean length of 16 bp. Mono-nucleotide repeats accounted for the largest percentage (77.05%) among all derived loci, followed by di-nucleotide (10.18%) and tetra-nucleotide (6.12%) repeats ( Table 1). The frequency of tri-nucleotide, penta-nucleotide and hexa-nucleotide repeats was lower, being 3.56, 2.14 and 0.95%, respectively. Our study found thousands of times more microsatellite sequences than those obtained by using the DNA library and the probe hybridization enrichment methods [31]. For example, only 58 microsatellite-containing fragments were obtained by screening genomic libraries with probes for B. alnoides [26] while 38 microsatellite sequences were obtained for B. pendula [32]. Among the abundant mono-to hexa-nucleotide repeats from the present study, mono-, di-and tetra-nucleotides were the major genomic microsatellite nucleotide types, while tri-, penta-and hexa-nucleotides were of relatively low proportion. Similar results have also been observed in other woody plants based on the same sequencing technology, such as Prunus persica, Carica papaya and Citrus sinensis [33].

Characteristics of SSR Sequences
It was found that 106,526 mono-nucleotide repeats (84.11% of the total) fell under 9-14 repeat number counts, and only four greater than 33 (Table 1). All the other nucleotide repeats were

Distribution of SSR in the Genome
A total of 64,376 contig sequences, about 60.47%, of 106,452 clean reads, contained SSR loci. Among these, 25,431 sequences contained only one SSR locus, and the other 38,945 sequences contained at least two SSR loci. In total 164,357 candidate SSR loci were obtained ( Table 1). The total size of examined sequences was 106,901,537 bp, and the size of the SSR repeat lengths ranged from 10 bp to 75 bp with a mean length of 16 bp. Mono-nucleotide repeats accounted for the largest percentage (77.05%) among all derived loci, followed by di-nucleotide (10.18%) and tetra-nucleotide (6.12%) repeats (Table 1). The frequency of tri-nucleotide, penta-nucleotide and hexa-nucleotide repeats was lower, being 3.56, 2.14 and 0.95%, respectively. Our study found thousands of times more microsatellite sequences than those obtained by using the DNA library and the probe hybridization enrichment methods [31]. For example, only 58 microsatellite-containing fragments were obtained by screening genomic libraries with probes for B. alnoides [26] while 38 microsatellite sequences were obtained for B. pendula [32]. Among the abundant mono-to hexa-nucleotide repeats from the present study, mono-, di-and tetra-nucleotides were the major genomic microsatellite nucleotide types, while tri-, penta-and hexa-nucleotides were of relatively low proportion. Similar results have also been observed in other woody plants based on the same sequencing technology, such as Prunus persica, Carica papaya and Citrus sinensis [33].

Characteristics of SSR Sequences
It was found that 106,526 mono-nucleotide repeats (84.11% of the total) fell under 9-14 repeat number counts, and only four greater than 33 (Table 1). All the other nucleotide repeats were mainly within 3-8 repeat number category, 82.37% for di-, 96.91% for tri-, and nearly 100% for tetra-, pentaand hexa-nucleotide repeats, respectively.
There were 291 types of repeat motifs in total. AG/CT (46.33%) and AT/AT (44.15%) were the most common repeat motifs in four types of di-nucleotide, whereas CG/CG accounted for only 0.27% (Figure 2a). These results are in accordance with previous studies reported by Gupta et al. [34] and Powell et al. [35]. In addition, the GC motif of lower frequency may be associated with cytosine methylation, which can be converted into a thymidine by deamination [36]. Among the ten types of repeat motifs in tri-nucleotide found in the present study (Figure 1b), the highest frequency repeat motif was AAT/ATT (48.98%), and the lowest was CCG/CGG (0.72%, Figure 2b). The main repeat motifs of tri-nucleotide were different from other plant species such as maize (CCG/GGC), rice (AGG/TCC) and barley (CCG/GGC) [37]. Repeat motifs AAAT/ATTT (60.85%), AAAAT/ATTTT (32.29%) and AAAAAT/ATTTTT (25.91%) observed in the present study were the most frequent in tetra-, pentaand hexa-nucleotide, respectively (Figure 2c-e).
Molecules 2018, 23, x FOR PEER REVIEW 4 of 10 mainly within 3-8 repeat number category, 82.37% for di-, 96.91% for tri-, and nearly 100% for tetra-, penta-and hexa-nucleotide repeats, respectively. There were 291 types of repeat motifs in total. AG/CT (46.33%) and AT/AT (44.15%) were the most common repeat motifs in four types of di-nucleotide, whereas CG/CG accounted for only 0.27% (Figure 2a). These results are in accordance with previous studies reported by Gupta et al. [34] and Powell et al. [35]. In addition, the GC motif of lower frequency may be associated with cytosine methylation, which can be converted into a thymidine by deamination [36]. Among the ten types of repeat motifs in tri-nucleotide found in the present study (Figure 1b), the highest frequency repeat motif was AAT/ATT (48.98%), and the lowest was CCG/CGG (0.72%, Figure 2b). The main repeat motifs of tri-nucleotide were different from other plant species such as maize (CCG/GGC), rice (AGG/TCC) and barley (CCG/GGC) [37]. Repeat motifs AAAT/ATTT (60.85%), AAAAT/ATTTT (32.29%) and AAAAAT/ATTTTT (25.91%) observed in the present study were the most frequent in tetra-, penta-and hexa-nucleotide, respectively (Figure 2c-e).

Primer Screening and SSR Markers Polymorphism Detecting
In total 2549 pairs of SSR primers were designed from all candidate SSR loci, of which 138 pairs were excluded due to the occurrence of compound formation in repeat regions of relevant SSR sequences. The majority of nucleotide repeats were di-nucleotide (75.53%), followed by tri-nucleotide (20.24%). A total of 900 pairs of primers were randomly selected for evaluating the successfulness of amplification, among them di-, tri-, tetra-, penta-and hexa-nucleotide repeats accounted for 74.11%, 21.00%, 3.33%, 0.67%, and 0.89%, respectively, and 580 pairs of primers (64.44%) were amplified successfully with clear and distinguishable microsatellite bands. The optimum annealing temperature of these SSR primers was 60 • C as determined by a temperature trial with three individuals. Eventually, 310 primers were shown to detect polymorphisms in 24 samples from a natural B. alnoides forest in Jingxi County, Guangxi, China and are therefore designated polymorphic SSR loci. Characterizations of 310 polymorphic SSR loci are presented in Table S1, and the remaining 270 pairs generated monomorphic markers or fixed heterozygosity in 24 samples or had weak fluorescence signals in amplified products of most samples (Table S2). The fragment length of the amplified products was from 105 bp to 306 bp, which was within the expectation range (100-400 bp). None of the 310 SSR markers developed in the present study was the same as the 19 pairs developed in a previous study [26]. When we tested the usefulness of the 19 pairs developed by Guo et al. [26] for our 24 B. alnoides samples, it was found that all of them were identified as polymorphic loci (Table S3). Based on the previously reported SSR markers (Figure 3a) and the SSR markers developed in the present study (Figure 3b), a phylogenetic relationship was constructed between the 24 individuals of B. alnoides. It could be seen that the larger number of markers resulted in very different phylogenetic relationships between the 24 samples.

Primer Screening and SSR Markers Polymorphism Detecting
In total 2549 pairs of SSR primers were designed from all candidate SSR loci, of which 138 pairs were excluded due to the occurrence of compound formation in repeat regions of relevant SSR sequences. The majority of nucleotide repeats were di-nucleotide (75.53%), followed by tri-nucleotide (20.24%). A total of 900 pairs of primers were randomly selected for evaluating the successfulness of amplification, among them di-, tri-, tetra-, penta-and hexa-nucleotide repeats accounted for 74.11%, 21.00%, 3.33%, 0.67%, and 0.89%, respectively, and 580 pairs of primers (64.44%) were amplified successfully with clear and distinguishable microsatellite bands. The optimum annealing temperature of these SSR primers was 60 °C as determined by a temperature trial with three individuals. Eventually, 310 primers were shown to detect polymorphisms in 24 samples from a natural B. alnoides forest in Jingxi County, Guangxi, China and are therefore designated polymorphic SSR loci. Characterizations of 310 polymorphic SSR loci are presented in Table S1, and the remaining 270 pairs generated monomorphic markers or fixed heterozygosity in 24 samples or had weak fluorescence signals in amplified products of most samples (Table S2). The fragment length of the amplified products was from 105 bp to 306 bp, which was within the expectation range (100-400 bp). None of the 310 SSR markers developed in the present study was the same as the 19 pairs developed in a previous study [26]. When we tested the usefulness of the 19 pairs developed by Guo et al. [26] for our 24 B. alnoides samples, it was found that all of them were identified as polymorphic loci (Table S3). Based on the previously reported SSR markers (Figure 3a) and the SSR markers developed in the present study (Figure 3b), a phylogenetic relationship was constructed between the 24 individuals of B. alnoides. It could be seen that the larger number of markers resulted in very different phylogenetic relationships between the 24 samples.    Table S1. The number of alleles (N A ) of each locus ranged from 2 to 19 with a mean of 5.14, and 86.45% of loci had 2 to 8 alleles. The observed (H O ) and expected (H E ) heterozygosities of each locus varied from 0.04 to 1.00 (mean 0.64) and 0.04 to 0.92 (mean 0.57) respectively. The Shannon-Wiener diversity index (I) ranged from 0.10 to 2.68 with a mean of 1.12. The parameters of polymorphism decreased as nucleotide repeat varying from di-to tetra-nucleotides ( Figure 4). However, this trend could not be seen for penta-and hexa-nucleotides due perhaps to their less occupancy. As a whole, the percentage of polymorphic loci with di-nucleotide repeats was larger than that with other nucleotide repeats. It was indicated that a selection of di-nucleotide repeats might be more effective to screen out polymorphic microsatellite markers, which is similar to the results of a previous study on B. platyphylla [38]. The genetic variation revealed by all polymorphic SSR loci for 24 individuals is also demonstrated in Table S1. The number of alleles (NA) of each locus ranged from 2 to 19 with a mean of 5.14, and 86.45% of loci had 2 to 8 alleles. The observed (HO) and expected (HE) heterozygosities of each locus varied from 0.04 to 1.00 (mean 0.64) and 0.04 to 0.92 (mean 0.57) respectively. The Shannon-Wiener diversity index (I) ranged from 0.10 to 2.68 with a mean of 1.12. The parameters of polymorphism decreased as nucleotide repeat varying from di-to tetra-nucleotides ( Figure 4). However, this trend could not be seen for penta-and hexa-nucleotides due perhaps to their less occupancy. As a whole, the percentage of polymorphic loci with di-nucleotide repeats was larger than that with other nucleotide repeats. It was indicated that a selection of di-nucleotide repeats might be more effective to screen out polymorphic microsatellite markers, which is similar to the results of a previous study on B. platyphylla [38].

Cross-Species Transferability
Ninety six pairs of SSR primers were randomly selected to examine cross-species transferability among six related species in the genus Betula. The results showed that 48.96% of these SSR primers could be amplified in B. platyphylla, 73.96% in B. austro-sinensis, 82.29% in B. cylindrostachya, 78.13% in B. fujianensis, 79.17% in B. hainanensis, and 84.38% in B. luminifera (Table  S4). These results thus indicated that these SSR primers had a high rate of transferability among the six species. High rate of transferability between species in the same genus has also been verified by Wang et al. [39] who discovered that 86.10% of SSR primer pairs developed from other Prunus species was transferable to P. virginiana. This is because the transferability of SSR primers across species is determined by the conservativeness of flanking sequences of microsatellite and their evolutionary stability [40]. The genomic differences between the related species are relatively small, making SSR primers highly transferable across specie. Thus, a large number of available primer NA HO HE I

Cross-Species Transferability
Ninety six pairs of SSR primers were randomly selected to examine cross-species transferability among six related species in the genus Betula. The results showed that 48.96% of these SSR primers could be amplified in B. platyphylla, 73.96% in B. austro-sinensis, 82.29% in B. cylindrostachya, 78.13% in B. fujianensis, 79.17% in B. hainanensis, and 84.38% in B. luminifera (Table S4). These results thus indicated that these SSR primers had a high rate of transferability among the six species. High rate of transferability between species in the same genus has also been verified by Wang et al. [39] who discovered that 86.10% of SSR primer pairs developed from other Prunus species was transferable to P. virginiana. This is because the transferability of SSR primers across species is determined by the conservativeness of flanking sequences of microsatellite and their evolutionary stability [40]. The genomic differences between the related species are relatively small, making SSR primers highly transferable across specie. Thus, a large number of available primer pairs can be used in the genetic studies and molecular assistant breeding for the species in the genus Betula.

Plant Material and DNA Isolation
In total 24 individuals of B. alnoides were sampled from natural forest in Jingxi County, Guangxi, China. Among them, one sample was randomly selected for de novo genome sequencing, three samples were applied to screen out the successfully amplified SSR primers and optimum annealing temperatures, and all 24 individuals were used to evaluate polymorphism of the developed SSR loci. Six species in the genus Betula (Table 2) were used to assess cross-species transferability of these SSR loci. Fresh leaves were collected from each individual and dried separately with silica gel. Total genomic DNA was then extracted by modified CTAB method [41] and stored at −20 • C. The concentration and purity of DNA were determined by NanoDrop 2000 (Thermo Fisher Scientific Inc., Waltham, MA, USA).

Genome Sequencing, SSR Finding and Survey
De novo genome sequencing for B. alnoides was conducted using 1/6 of a run on the Roche 454 GS FLX + platform (454 Life Sciences, Roche Company, Branford, CT, USA) according to the process of Guo et al. [42]. For each sequence of raw data, the quality-control processes were conducted with the software Qiime (version 1.17, http://qiime.org/), which included: (1) the sequencing adapters were trimmed; (2) low-quality bases were trimmed with the index of the Q20 bases percent more than 90%; and (3) the ambiguous bases were deleted.
The MISA software was used to search SSR with following conditions: (1) the repeat number of the mono-nucleotide was at least ten; (2) the repeat number of the di-nucleotide was at least five; (3) the repeat number of the tri-nucleotide was at least four; and (4) the repeat numbers of the te-, penta-, and hexa-nucleotide were at least three. The characteristics of SSR were analyzed at the whole genome level.
The optimum annealing temperature of SSR primers was screened out from 56 • C, 58 • C and 60 • C. The PCR reaction mixture (10 µL) contained 50 ng of DNA template, 150 µM dNTPs, 2.0 µM MgCl 2 , 0.5 µM forward and reverse primers, 1× PCR buffer (Tiangen Biotech Ltd., Beijing, China), and 0.04 U/µL of Taq DNA polymerase (Tiangen Biotech Ltd., Beijing, China). PCR was performed on the PCR system (Applied Biosystems Veriti) according to the following program: 4 min for initial denaturation at 94 • C; 31 cycles of 30 s for denaturation at 94 • C, 30 s at annealing temperature and 30 s at 72 • C; and 10 min at 72 • C. The amplified products were stored at 4 • C, and detected by 1% agarose gel electrophoresis. Successful amplification was judged by clear and distinguishable microsatellite bands.

SSR Markers Polymorphism Detection
Polymorphism of SSR loci was detected with the M13-tailed primers method [43]. SSR forward primers were synthesized and labeled with M13 sequence (5 -CACGACGTTGTAAAACGAC-3 ) at 5 end, and the reverse primers were not labeled. The M13 primer was labeled with a fluorescent dye FAM, NED, VIC, or ROX. The PCR was performed with a reaction mixture (10 µL) containing 50 ng of DNA template, 150 µM dNTPs, 2.0 µM MgCl 2 , 0.5 µM fluorescent forward M13-labeled primers and reverse primers, 1× PCR buffer (Tiangen Biotech Ltd., Beijing, China), and 0.04 U/µL of Taq DNA polymerase (Tiangen Biotech Ltd., Beijing, China). The PCR program was used as mentioned above.
The PCR products were detected using 3730 XL automatic sequencer (ABI Co., Foster, CA, USA). The genotyping of products was performed using GeneMarker V2.2.0 [44]. The number of alleles (N A ) and observed heterozygosity (H O ) were calculated with the results of genotyping in Excel 2010 (Microsoft Corporation, Redmond, WA, USA), the expected heterozygosity (H E ) and Shannon-Wiener diversity index (I) were generated using ATETRA [45], which is specially applied for the analysis of tetraploid. The Nei's genetic distance between the 24 samples was computed by NTSYS 2.1 [46] software. The software was also used to obtain clustering graph of 24 samples based on Nei's genetic distance by the unweighted pair group method with arithmetic (UPGMA) method.

Conflicts of Interest:
The authors declare no conflict of interest.