extfieldsofstudy: Medicine Biology

text:
Role of Antibodies in Protection Elicited by Active Vaccination with Genetically Inactivated Alpha Hemolysin in a Mouse Model of Staphylococcus aureus Skin and Soft Tissue Infections

ABSTRACT Due to the emergence of highly virulent community-associated methicillin-resistant Staphylococcus aureus (CA-MRSA) infections, S. aureus has become a major threat to public health. A majority of CA-MRSA skin and soft tissue infections in the United States are caused by S. aureus USA300 strains that are known to produce high levels of alpha hemolysin (Hla). Therefore, vaccines that contain inactivated forms of this toxin are currently being developed. In this study, we sought to determine the immune mechanisms of protection for this antigen using a vaccine composed of a genetically inactivated form of Hla (HlaH35L). Using a murine model of skin and soft tissue infections (SSTI), we found that BALB/c mice were protected by vaccination with HlaH35L; however, Jh mice, which are deficient in mature B lymphocytes and lack IgM and IgG in their serum, were not protected. Passive immunization with anti-HlaH35L antibodies conferred protection against bacterial colonization. Moreover, we found a positive correlation between the total antibody concentration induced by active vaccination and reduced bacterial levels. Animals that developed detectable neutralizing antibody titers after active vaccination were significantly protected from infection. These data demonstrate that antibodies to Hla represent the major mechanism of protection afforded by active vaccination with inactivated Hla in this murine model of SSTI, and in this disease model, antibody levels correlate with protection. These results provide important information for the future development and evaluation of S. aureus vaccines.

S taphylococcus aureus is an opportunistic pathogen that colonizes 30% of the population (1). This bacterium has a plethora of virulence factors that it uses to evade the immune system, scavenge for needed nutrients, and attach to and destroy host tissue (2,3). S. aureus is a common cause of uncomplicated skin infections (4) but can also cause more invasive infections, like sepsis, endocarditis, and osteomyelitis (2). Methicillin-resistant S. aureus (MRSA), which is resistant to all ␤-lactam antibiotics, has historically been associated with nosocomial infections; however, community-acquired MRSA (CA-MRSA) has emerged in recent years (5). CA-MRSA infections, a majority of which in the United States are caused by the highly virulent USA300 genotype (6,7), often occur in immunocompetent patients lacking obvious risk factors (4). Vancomycin is the drug of choice for treating MRSA infections (8,9). However, due to its extensive use (5,9), vancomycinresistant S. aureus (VRSA) has been increasingly isolated since the first reported vancomycin-intermediate S. aureus (VISA) infection in 1997 (10,11). The emergence of hypervirulent CA-MRSA and VRSA strains is an alarming threat to public health that highlights the need for an S. aureus vaccine.
Others have suggested that prophylactic therapies that target CA-MRSA should include inactivated alpha hemolysin (Hla), since Hla is believed to contribute to the severity of USA300 skin infections (12) and is expressed at high levels by USA300 strains (13). Hla is a secreted toxin that binds to the eukaryotic cell surface protein ADAM10 and then oligomerizes to generate a heptameric pore that is cytolytic to multiple cell types, including erythrocytes and endothelial cells (14). Passive immunization with anti-Hla monoclonal antibodies has been shown to reduce lesion size in a mouse dermonecrosis model (15). Inactivated forms of Hla have been included in vaccine candidates that have been studied in the clinic setting (16). A better understanding of the mechanism(s) of protection for this antigen would improve our ability to evaluate these vaccines. For example, information re-garding the immune mechanisms of protection for Hla-based vaccines ensures that appropriate immune responses to vaccination are measured both in preclinical and clinical studies. Moreover, the development of correlates of protection can be useful for monitoring vaccine consistency and the susceptibilities of specific populations postvaccination (17).
In this study, we examined the mechanisms of protection of a vaccine based on the genetically inactivated alpha hemolysin mutant HlaH35L, in which histidine 35 was substituted with leucine, in a mouse model of S. aureus skin and soft tissue infections (SSTI). HlaH35L binds to the cell membrane; however, it is unable to form a pore and is therefore not toxic (18,19). We determined the important role that antibodies play in the protection afforded by active immunization with this genetically inactivated Hla-based vaccine and demonstrate that anti-Hla antibody levels correlate with protection.

Ethics statement.
For all mouse studies, the protocols were reviewed and approved by the Institutional Animal Care and Use Committee (IACUC) of the Center for Biologics Evaluation and Research.

Reagents.
Purified recombinant HlaH35L was generated as described by Brady et al. (20).
Rabbit anti-HlaH35L polyclonal antibody production and IgG purification. The production of Rabbit anti-HlaH35L polyclonal antibody was performed by Cocalico Biologicals, Inc. (Reamstown, PA), subject to their IACUC approval. New Zealand White rabbits (2 to 3 months old, male, 1.8 to 2.3 kg) were immunized with 100 g HlaH35L at day 0. The rabbits were then immunized again on days 21 and 42 with 50 g HlaH35L. For each dose, the antigen was diluted to 0.5 ml in 0.9% saline and adsorbed to an equal volume of Alhydrogel (Brenntag Biosector, Frederikssund, Denmark). Bleeds were obtained from the rabbits at days 56 and 70. IgG was purified by affinity chromatography. Rabbit preimmune and immune sera were diluted 1:1 in protein A IgG binding buffer (provided by Thermo Scientific, Rockford, IL) and applied to 5 ml protein A columns (NAb Spin columns, Thermo Scientific, Rockford, IL). The columns were placed on a rotator and incubated for 10 min at room temperature. The columns were then washed three times with 10 ml of binding buffer, and IgG was eluted three times with 5 ml of IgG elution buffer (provided by Thermo Scientific, Rockford, IL), followed by immediate neutralization of the pH using 0.5 ml of 1 M Tris-HCl (pH 8) (Fisher Scientific, Fair Lawn, NJ). Pooled fractions were diafiltered into phosphate-buffed saline (PBS) (Quality Biological, Gaithersburg, MD) using Amicon Ultra-15-ml 30,000 molecular weight cutoff (MWCO) centrifugal filters (EMD Millipore, Billerica, MA). IgG concentration was measured using the GE NanoVue (GE Healthcare, Piscataway, NJ) IgG function (which measures the A 280 and calculates protein concentration using the extinction coefficient of IgG). The purified IgG was characterized on NuPAGE 4 to 12% Bis-Tris gels in morpholinepropanesulfonic acid (MOPS) buffer (Invitrogen, Grand Island, NY). The final IgG concentration was brought to 22 to 24 mg/ml. Active immunization studies. Purified recombinant HlaH35L was thawed and diafiltered into PBS using Amicon Ultra-15-ml 30,000 MWCO centrifugal filters. The HlaH35L concentration was determined by the bicinchoninic acid (BCA) Protein Assay (Pierce, Rockford, IL) according to the manufacturer's protocol. HlaH35L was then adsorbed to Alhydrogel with CpG for 1 h on ice (CBER Core Facility, Bethesda, MD) to obtain final concentrations of 200 g/ml HlaH35L, 2 mg/ml Alhydrogel, and 150 g/ml CpG. We used this antigen/adjuvant combination because we have shown that this formulation is protective in a murine SSTI model of S. aureus infection. Six-week-old female BALB/c mice or Jh mice were immunized subcutaneously in the flank with 100 l formulated vaccine containing 20 g of adjuvanted HlaH35L on days 0 and 14. All mice were purchased from Taconic, Hudson, NY. The control mice received 100 l of adjuvant only, which was prepared as described above but without HlaH35L. The mice were bled at day 26 to determine antibody levels and were challenged with S. aureus at day 28, as described below. The dose-response study was conducted as described above except that 6-week-old female BALB/c mice (NCI, Frederick, MD) were immunized only once (day 0), and the highest vaccine dose was prepared as described above at a volume sufficient to allow for 1:4 serial dilutions following adsorption. Extra control (adjuvant only) vaccine was prepared and used as a diluent in these serial dilutions.
Passive immunization study. Eight-week-old female BALB/c mice (NCI, Frederick, MD) were passively immunized intraperitoneally with 200 l of concentrated IgG at days 0 and 3. The mice were challenged on day 1 with S. aureus, as described below.
Skin and soft tissue challenge. Two days prior to challenge, S. aureus strain SAP149 (21) was streaked on tryptic soy agar (TSA) with 5% sheep blood (BD, Franklin Lakes, NJ). The following day, an overnight culture was initiated using a single hemolytic colony. On the day of challenge, four baffled Erlenmeyer flasks (Corning, Corning, NY) containing 50 ml of tryptic soy broth (TSB) (MP Biomedicals, Solon, OH) were inoculated at a dilution of 1:50 using the overnight growth. The cultures were then placed in a 37°C shaking water bath and incubated until they reached an absorbance at 600 nm (A 600 ) between 0.75 and 0.8. The bacteria were pelleted by centrifuging at 4,000 ϫ g for 10 min. All four pellets were resuspended in 50 ml of PBS and the concentration was determined using a Petroff-Hausser counting chamber (Hausser Scientific Partnership, Horsham, PA). The bacteria were centrifuged again and resuspended with PBS to a final concentration of 1 ϫ 10 11 CFU/ml. The concentration was verified by serial dilution and plating on TSA. The mice were anesthetized and challenged as previously described by Brady et al. (20). Seven days postchallenge, each mouse was euthanized, and its affected ears were imaged when indicated. The affected ears were then gently cleaned with alcohol pads (PDI, Orangeburg, NY), removed, and homogenized in sterile PBS. The CFU count per ear was determined by serial dilutions and plating on TSA. In addition, the two individuals scored the ears for tissue damage based on a scale of 1 to 3, and these scores were averaged. Mice that received a score of 1 had no missing ear pieces, no scabs, little to no redness, and no or small holes. Mice that received a score of 2 had no missing ear pieces, small to medium scabs, moderate redness, and small to medium holes. Mice that received a score of 3 had missing pieces or large holes, medium to large scabs, and moderate to severe redness.
Enzyme linked immunosorbent assay. Enzyme-linked immunosorbent assays (ELISAs) to measure mouse antibodies were performed as described by Brady et al. (20). ELISAs to measure rabbit antibodies were performed according to the same protocol; however, the conjugate antibody used was goat anti-rabbit IgG conjugated to horseradish peroxidase (HRP) (Bio-Rad Laboratories, Hercules, CA). To determine IgG titers, the data were plotted using GraphPad Prism software (version 5, Graph-Pad Software, La Jolla, CA) with a 4-parameter logistic (4PL) curve fit of A 405 (y axis) versus the logarithm of the reciprocal serum dilution (x axis). The ELISA titers were defined as the 50% effective dose (ED 50 ) of a 4PL curve.
Hla neutralization assay. The assay was designed based on one described by Adhikari et al. (22). Briefly, serial dilutions of mouse serum were preincubated with 20 ng Hla (List Biological Laboratories, Campbell, CA) for 30 min at 37°C in 100 l total volume of PBS supplemented with 2.5% fetal bovine serum (Gibco, Life Technologies, Grand Island, NY) in a 96-well plate (MaxiSorp; Thermo Fisher, Waltham, MA). An equal volume of 1% rabbit erythrocytes (CBER Animal Facility, Bethesda, MD) was then added to each well, and incubation at 37°C was continued for another 30 min. Each row of diluted antiserum included a no-toxin control that served as a blank, and every plate included a toxin-only control for lysis, as well as serially diluted rabbit anti-Hla polyclonal antibody (the same as used for passive studies; see above) as a positive control for neutralization. After incubation, the plates were centrifuged at 500 ϫ g for 5 min, and 100 l of the supernatant from each well was transferred to a new plate. The degree of erythrocyte lysis was determined by reading the A 416 on a VersaMax microplate reader (Molecular Devices, Sunnyvale, CA). For each antibody dilution, the corresponding blank reading was subtracted to remove background absorbance, and the corrected absorbance readings were plotted against the log of the reciprocal antibody dilution using a 4PL-curve fit determined using GraphPad Prism. The 50% neutralization titer (NT 50 ) is defined as the reciprocal dilution of serum that neutralizes 50% of the hemolytic activity of Hla in the assay.
Statistical analyses. GraphPad Prism was used for statistical analyses. To assess the correlation between antibody titers and protection, the ELISA ED 50 was plotted as a function of the CFU/ear, and a linear best-fit curve was used to determine the r 2 value and whether the slope of the line was significantly different than zero.

RESULTS
Passive immunization with rabbit anti-HlaH35L IgG decreases bacterial colonization. Previously, others demonstrated in a murine skin infection model that anti-Hla polyclonal antibodies can passively protect against dermonecrotic lesions (12). However, because the model used in that study involved subcutaneous administration of S. aureus rather than the epicutaneous administration method used in our model, we sought to confirm that antibodies can passively protect in the epicutaneous model of skin infections used in this study. The mice were passively immunized with purified rabbit anti-HlaH35L IgG and normal rabbit IgG and then challenged using the epicutaneous model of skin infections developed by Prabhakara et al. (23). We found that passive immunization with anti-Hla polyclonal antibodies protected mice against lesions (data not shown), as was reported for the subcutaneous model (12). Moreover, passively immunized mice had significantly lower CFU counts/ear than control mice that received normal rabbit IgG (Fig. 1), indicating that passive immunization with anti-HlaH35L IgG not only protects against tissue damage but also reduces bacterial loads. These data indicate that antibodies are sufficient for protection in an epicutaneous model of S. aureus disease. However, the data do not exclude the possibility that other immune responses might contribute to protection or even be the overriding mechanisms of protection in an active immunization model. B-cell immunity is necessary in order to decrease disease severity and bacterial colonization. We therefore sought to determine whether antibodies are necessary for protection mediated by active immunization with an inactivated Hla-based vaccine, or whether other immune responses, such as T-cell responses, might represent the major protective mechanism(s). To this end, we examined the ability of genetically inactivated Hla to protect Jh mice, which are devoid of mature B lymphocytes and have no detectable IgM or IgG in their blood serum (24) against S. aureus infection in the SSTI model. Jh mice were immunized with either HlaH35L formulated with adjuvant or adjuvant only and were subsequently challenged with S. aureus. The effect of vaccination on the protection of Jh mice compared to the protection of wildtype BALB/c mice was then assessed. As was observed previously (20), at 14 days postchallenge, BALB/c mice immunized with adjuvanted HlaH35L had significantly less bacterial colonization ( Fig. 2A) and visibly less tissue damage of the affected ear (Fig. 2B) than did BALB/c mice that received the adjuvant control. In contrast, HlaH35L-immunized Jh mice did not exhibit a significant difference in bacterial colonization or in damage to the affected ear compared to control Jh mice immunized with adjuvant only ( Fig. 2A and B). These results demonstrate that B-cell-based immunity is necessary in order to decrease disease severity and bacterial colonization in the SSTI model.
Anti-HlaH35L antibody levels correlate with reduction in bacterial colonization. In order to determine if antibody levels correlate with bacterial colonization in the SSTI model, the mice were immunized once with either various doses of HlaH35L vaccine or adjuvant control, bled 2 days before challenge to measure antibody levels, and then were challenged with S. aureus. A clear dose response was observed, with significant protection achieved with a dose of HlaH35L vaccine as low as 0.3125 g (Fig. 3).
We next examined the correlation between antibody titer and bacterial levels at 14 days postchallenge. We observed a moderate correlation (r 2 ϭ 0.62) between total IgG titer (ED 50 ) and bacterial levels at 14 days postchallenge (Fig. 4A). The slope of the best-fit  line was significantly different from zero. These results indicate that anti-HlaH35L antibodies play a role in reducing bacterial colonization in the SSTI model, with a positive relationship between total IgG level and bacterial reduction. However, the moderate correlation suggests that total IgG titer may not fully explain the protective mechanisms of this antigen. Therefore, we also examined the relationship between toxin neutralizing antibody titer and reduction in bacterial levels. To do this, we examined the ability of serum samples from individual mice to neutralize Hla using a rabbit erythrocyte toxin neutraliza-tion assay. When these serum samples were titrated with fixed concentrations of erythrocytes and Hla, we noted a differential ability of individual samples to neutralize erythrocyte lysis. Unfortunately, due to the lack of sensitivity of the assay, we were unable to achieve full neutralization curves for all sera and were therefore not able to quantitatively determine the NT 50 values for a number of the serum samples. Instead, we assigned the neutralization capabilities of each the serum samples to one of three categories based on the extent of neutralization observed: (i) fully neutralizing (serum can completely inhibit hemolysis such that the neutralization curve reaches 0), (ii) partially neutralizing (the lowest absorbance reading of the neutralization curve is at less than one-half of the upper asymptote), and (iii) nonneutralizing (the lowest absorbance of the neutralization curve is more than one-half of the upper asymptote). We first qualitatively examined the relationship between vaccine dose and neutralization. As shown in Fig. 3, we observed that most animals that received higher doses of HlaH35L had anti-Hla antibodies that were at least partially neutralizing. We then examined which points on the ELISA ED 50 /CFU correlation best-fit curve fell into each neutralization category. The sera with neutralization capabilities (either fully neutralizing or partially neutralizing) tended to be correlated with higher ELISA titers and lower CFU levels in the SSTI model ( Fig. 4A; see also Table S1 in the supplemental material). As shown in Fig. 4B, mice that had either partially or fully neutralizing titers had significantly lower bacterial levels than control mice that received adjuvant alone. These results suggest that toxin neutralization may be important in protecting mice from SSTI infection after active immunization with HlaH35L.

DISCUSSION
Hla is a well-studied toxin, which, in its inactivated form, has been included in vaccines under development (16); however, the immune mechanisms of Hla-based vaccine-induced protection and immune parameters that correlate with protection for active vaccination with this antigen have not been precisely defined. In this study, we used a mouse model of S. aureus SSTI to examine the The results are representative of two independent experiments. The black symbols represent mice with nonneutralizing sera, blue symbols represent mice with partially neutralizing sera, and red symbols represent mice with fully neutralizing sera, as determined using the Hla neutralization assay described in Materials and Methods.  Fig. 3A. The slope of the line was determined to be significantly different from 0 (P Ͻ 0.0001), and the coefficient of determination (r 2 ) was calculated to be 0.62. The neutralizing titers were also determined for sera from each of the individual mice. The black symbols represent mice with nonneutralizing sera, blue symbols represent mice with partially neutralizing sera, and red symbols represent mice with fully neutralizing sera, as determined using the Hla neutralization assay described in Materials and Methods. (B) CFU counts/ear were compared between mice with partially neutralizing (), fully neutralizing (OE), and adjuvant control serum (). *, groups with a statistically significant decrease in CFU/ear compared to the adjuvant control group, as determined by 1-way ANOVA with a Tukey's multiple comparisons posttest; a P value of Ͻ0.05 was considered significant. The dashed line indicates the limit of detection for the assay. Each bar represents the mean Ϯ SEM. immune mechanisms of protection for a genetically inactivated Hla vaccine. We found that antibodies are both necessary and sufficient for protection in an epicutaneous SSTI model of infection. Antibodies were not only important for preventing tissue damage but also for reducing bacterial levels.
We further examined the relationship between antibodies and protection in order to determine whether antibodies levels correlate with protection. We found that IgG titers showed an inverse correlation with CFU levels; that is, higher antibody levels corresponded to decreased bacterial loads in affected ears. While the slope of the best-fit line was significantly nonzero, the moderate r 2 value of 0.62 suggests that total IgG levels do not fully explain the observed protection. One possibility for this observation is that only a subpopulation of antibodies confers protection. To that end, we attempted to correlate neutralizing antibody levels with protection; however, we were hampered in this effort because of the relative lack of sensitivity of the assays used to assess toxin neutralization. However, we did find that mice that had detectable neutralizing antibody levels were significantly protected, leaving open the possibility that very small amounts of neutralizing antibody might be sufficient to protect in this model. Such functional antibodies might represent the major mechanism of protection. Alternatively, other immune responses might contribute to protection; however, the major protective immune mechanism is clearly B-cell mediated, since Jh mice that are devoid of mature B lymphocytes but which have normal T-cell development (24) were not protected by active vaccination with a genetically inactivated Hla-based vaccine in this model.
Vaccine-induced protection is antibody mediated for a number of toxoid-based vaccines (17). Our work demonstrates that the protection that we observed using a genetically inactivated Hla vaccine in a murine SSTI model is likely similarly based. While we demonstrated that antibodies mediate protection in this SSTI model, the specific mechanism of how antibodies against a secreted toxin target S. aureus for pathogen clearance remains unknown. While a small amount of Hla is associated with the bacterial surface as it is secreted, the majority of toxin is not associated with the surface and, as such, it is unlikely that anti-Hla antibodies would opsonize the bacteria. Because neutrophils, which have been shown to be important in the control of S. aureus infections (23,(25)(26)(27)(28), are a target of alpha hemolysin (14), one might theorize that immunization with HlaH35L protects neutrophils, or other cells involved in the immune response that are important for control of the bacteria, from intoxication with Hla during infection. Alternatively, because secreted toxins cause localized tissue damage, the neutralization of alpha hemolysin may decrease subsequent nutrient release from host tissues and lessen the ability of S. aureus to further invade the host. This idea is supported by decreased levels of tissue destruction in HlaH35L-immunized mice (20). Whatever the mechanism, active toxin appears to be involved in the survival of the bacteria (23), and anti-Hla antibodies clearly play a pivotal role in the protection afforded by active immunization with an inactivated Hla-based vaccine in this SSTI model.
In summary, we determined that antibodies mediate protection afforded by vaccination with a genetically inactivated Hlabased vaccine in a murine model of S. aureus SSTI. Because inactivated Hla-based vaccines are currently undergoing clinical study, knowledge of the immune parameters that correlate with protection for such vaccines would ensure that the appropriate immune markers are evaluated in clinical trials. We believe that the work described here represents an important first step in the development of correlates of protection for new S. aureus vaccines that contain an inactivated Hla component. The work reported here was conducted using a murine model of S aureus SSTI; however, we recognize that what is true for mice might not necessarily be true for humans. In this regard, Fritz et al. (29) recently found that elevated titers to Hla were significantly correlated with protection from recurrent S. aureus infection in humans. These results are consistent with the findings reported here and underscore the importance of anti-Hla antibodies in protection. Additional clinical studies are needed to confirm the role that anti-Hla antibodies play in the protection induced by active immunization against S. aureus infections in humans and to establish a protective level.


extfieldsofstudy: Biology Medicine

text:
The Li2 Mutation Results in Reduced Subgenome Expression Bias in Elongating Fibers of Allotetraploid Cotton (Gossypium hirsutum L.)

Next generation sequencing (RNA-seq) technology was used to evaluate the effects of the Ligon lintless-2 (Li2) short fiber mutation on transcriptomes of both subgenomes of allotetraploid cotton (Gossypium hirsutum L.) as compared to its near-isogenic wild type. Sequencing was performed on 4 libraries from developing fibers of Li2 mutant and wild type near-isogenic lines at the peak of elongation followed by mapping and PolyCat categorization of RNA-seq data to the reference D5 genome (G. raimondii) for homeologous gene expression analysis. The majority of homeologous genes, 83.6% according to the reference genome, were expressed during fiber elongation. Our results revealed: 1) approximately two times more genes were induced in the AT subgenome comparing to the DT subgenome in wild type and mutant fiber; 2) the subgenome expression bias was significantly reduced in the Li2 fiber transcriptome; 3) Li2 had a significantly greater effect on the DT than on the AT subgenome. Transcriptional regulators and cell wall homeologous genes significantly affected by the Li2 mutation were reviewed in detail. This is the first report to explore the effects of a single mutation on homeologous gene expression in allotetraploid cotton. These results provide deeper insights into the evolution of allotetraploid cotton gene expression and cotton fiber development.

Introduction
Cotton is the major source of natural fibers used in the textile industry. There are four cultivated species: AA genome diploids, Gossypium arboretum L. and G. herbaceum L.; and AADD genome allotetraploids, G. hirsutum L. and G. barbadense L. Upland cotton (G. hirsutum) represents about 95% of world cotton production [1]. Allotetraploid species originated around 1-2 million years ago from inter-specific hybridization between an AA-genome diploid native to Africa and Mexican DD-genome diploid [1,2].
Cotton fibers are single-celled trichomes that emerge from the ovule epidermal cells. About 25-30% of the seed epidermal cells differentiate into spinnable fibers [3]. Fiber length ranges from short (fuzz ,6 mm) to long (lint). Lint fibers of economically important G. hirsutum generally grow up to about 30-40 mm in length. Cotton fiber development undergoes four distinctive but overlapping stages: initiation, elongation, secondary cell wall biosynthesis, and maturation [4]. The rate and duration of each developmental stage is important to the quality attributes of the mature fiber. Cell elongation is crucial for fiber length, whereas secondary cell wall thickening is important for fiber fineness and strength.
Cotton fiber mutants are useful tools to elucidate biological processes of cotton fiber development. A cotton plant with abnormally short lint fibers was discovered in a breeding nursery of the Texas Agricultural Experiment Station in 1984. This mutant had short lint fibers (,6 mm) visually similar to those produced by Ligon lintless-1 (Li 1 ); however, unlike the stunted and deformed vegetative morphology caused by the Li 1 mutation, this fiber mutant had normal vegetative growth. The trait was controlled by one dominant gene named Ligon lintless-2 (Li 2 ) [5]. This gene was mapped to chromosome 18 (D T subgenome of G. hirsutum) using several approaches [6][7][8]. In a fiber developmental study, Kohel et al. [9] observed that elongation is restricted in Li 2 fibers, however secondary wall development proceeds normally in proportion to fiber length. Two near-isogenic lines (NILs) of Li 2 with the Upland cotton variety DP5690 were developed in a backcross program at Stoneville, MS [6]. Morphological evaluation of developing fibers did not reveal apparent differences between WT and Li 2 NILs during initiation or early elongation up to 5 days post-anthesis (DPA). Transcript and metabolite evaluations revealed significant changes in biological processes associated with cell expansion in the Li 2 mutant line at peak of fiber elongation, including reactive oxygen species, hormone homeostasis, nitrogen metabolism, carbohydrate biosynthesis, cell wall biogenesis, and cytoskeleton [6,10]. Therefore, the Li 2 mutation can be considered as a factor affecting cotton fiber elongation process, making it an excellent model system to study cotton fiber elongation.
In previous reports, we used microarray techniques to investigate global gene expression in Li 2 NILs [6,10]. However, by using the genome sequence of G. raimondii [11], RNA-seq can provide a more comprehensive and accurate transcriptome analysis based on the reference DNA sequences [12]. RNA-seq offers a larger dynamic range of quantification, reduced technical variability, and higher accuracy for distinguishing and quantifying expression levels of homeologous copies than DNA microarrays [12]. Because of the limited sequence divergence between the A T and D T subgenomes in cotton [13], a pipeline was developed to map and categorize RNA-seq reads as originating from the A T or D T subgenomes [14].
In the present study we compared quantitative gene expression levels of RNA-seq data between developing fibers of Li 2 and its WT NILs. We investigated the Li 2 mutation's effect on global transcriptional changes in subgenomes and on the functional distribution of homeologous genes during fiber elongation. These results provide deeper insights into the evolution of allotetraploid cotton gene expression.

RNA-seq of Wild Type and Li 2 Developing Fibers at Peak of Elongation
Considering the cost of deep sequencing, only one time point, at the peak of elongation, was selected for RNA-seq analysis, including two biological replicates for wild type and mutant NILs. The time points 8 and 12 days post-anthesis (DPA) represent peak rates of fiber elongation. The time point 8 DPA was selected because: 1) our earlier research revealed significant transcript and metabolite changes between the Li 2 and wild type NILs during this time of fiber development [6,10]; 2) the transcript level of the elongation stage-related gene GhExp1 significantly decreased in Li 2 mutant fiber at 8 DPA [6].
A total of 639 million reads (each 101 bp in length) from 4 libraries were obtained by paired-end Illumina sequencing. Approximately 2.3% more reads were obtained from Li 2 than wild type fiber transcriptomes. From 84.4% to 90.2% of the reads were mapped to the D 5 -genome reference sequence of G. raimondii (Table 1). Not all the reads mapped to the reference genome sequence, probably since some of the genes were not included in the 13 large pseudo-molecules and transcripts mapped to genomic regions were outside of annotated genes. Of the mapped reads, between 29.3%-31.4% were mapped to the A T subgenome and between 23.4%-25.1% were mapped to the D T subgenomes of G. hirsutum. If the mapped reads overlapped a homeologous SNP position (SNPs between the A T and D T subgenomes), they were categorized as belonging to one of the two subgenomes or as a chimeric read (A-reads, D-reads, and X-reads, respectively; [14]). If a read did not overlap a homeologous SNP position, the read was unable to be categorized as originating from either the A T or D T subgenome (N reads; Table 1). Notably, more reads from each library were aligned to the A T subgenome than to the D T subgenome. Among the 37,223 genes on the 13 chromosomes of the G. raimondii genome, 34,692 genes (93%) had at least one mapped read from developing fibers at peak of elongation (Table 2).

Differential Gene Expression in Developing Fibers
Counts of mapped reads were evaluated in wild type and mutant fiber transcriptomes. Genes were considered to be expressed if they had $10 reads mapped in one sample. Genes that were not considered to be expressed were not included in further analyses. Approximately 3% more expressed genes were detected in Li 2 than in wild type. Of the 37,223 genes on the 13 chromosomes of the G. raimondii D 5 reference genome 29,603 (79.5%) genes were expressed in wild type and 30,842 (82.9%) genes were expressed in Li 2 fiber (Table 2).
Many genes had altered expression levels as a result of the Li 2 mutation (Table S1in File S1). Some genes were expressed in one treatment (such as Li 2 ) but not the other treatment. For example, expressions of genes annotated as SAUR-like auxin-responsive protein (Gorai.005G257000), bHLH (Gorai.003G034700) and NAC domain transcription factor (Gorai.009G170700) were only detected in wild-type fiber. Cytokinin response factor 6 (Gorai.007G105600), UGT73C14 (Gorai.002G107900), cystein proteinase (Gorai.007G329600), MYB-like 102 (Gorai.012G132200) and WRKY transcription factor (Gorai.009G157300) were only detected in Li 2 fiber. The majority of these genes have not yet been functionally characterized in cotton, except of glycosyltransferase UGT73C14, which has been shown recently to be involved in ABA homeostasis [15].
The quantitative levels of mapped reads were evaluated for differential expression between elongating fibers of Li 2 and wild type NILs. A gene-by-gene ANOVA determined that 7,163 of 31,114 expressed genes were differentially expressed (FDR corrected p-value ,0.05) and had $2-fold difference in at least one of the following comparisons: fiber type (Li 2 versus wild type), A T /D T subgenomes, and combinations of these factors (statistical data for significantly regulated genes are provided in Data S1). The highest numbers of significantly differentially expressed genes were identified between the A T and D T subgenomes in wild type and Li 2 ; whereas approximately 3 times fewer differentially expressed genes were detected between fiber type comparisons ( Figure 1A). Of the 29,603 expressed homeologous pairs in wild type, 4,578 (wtA/wtD, 15.5%) showed significantly different expression level between subgenomes; whereas in mutant fiber of the 30,842 expressed genes, 3,967 (LiA/LiD, 12.9%) were differentially expressed between subgenomes. Therefore, the homeolog expression bias was significantly (p-value ,0.0001; Chi square) reduced in Li 2 fiber transcriptome.
In general, the A T subgenome contributed more differentially expressed genes to fiber transcriptome than did the D T subgenome. Approximately two times more genes were differentially expressed in the A T subgenome compared to the D T subgenome in wild type (A T -2958 vs. D T -1620; Figure 1B) and mutant fiber (A T -2574 vs. D T -1393). Comparison between fiber types showed more genes were upregulated in Li 2 versus wild type in both subgenomes ( Figure 1B). It should be noted that only about 38% (583 genes out of 1,536 in A T and 1,511 in D T subgenome) of significantly regulated genes between mutant and wild type overlapped between subgenomes ( Figure 1A).

Mutation Effects on Transcriptome of A T and D T Subgenomes of Allotetraploid G. hirsutum
The effect of Li 2 mutation on the transcriptome of each subgenome was evaluated. The genes significantly (FDR corrected p-value ,0.05) up-regulated ($2-fold) in one subgenome of wild type were considered to have biased expression. Of the 2,958 A T biased genes, 26.5% (784) had significantly changed the expression levels in both subgenomes of Li 2 as a result of mutation. However, of the 1,620 D T biased genes, 35.9% (582) had significantly changed the expression levels in both subgenomes of the Li 2 mutant ( Figure 2). Therefore Li 2 has a significantly greater effect (p-value ,0.0001; Chi square) on D T biased genes than A T biased genes. Importantly, the majority of biased genes had significantly reduced expression levels (8.6% and 12.4%), whereas only a small portion of genes increased expression levels (1.6% and 2.7%) in the mutant. However, more genes, which were down-regulated in homeologous subgenome in wild type, had increased expression levels in the mutant fiber: 11.9% and 14% were up-regulated, whereas 4.4% and 6.8% were down-regulated ( Figure 2).
Furthermore, a few homeolog pairs had reciprocal expression biases between two subgenomes as a result of mutation. Expression levels for three of these genes were tested by RT-qPCR across eight developmental time points from DOA to 20 DPA, representing initiation, elongation, and beginning of secondary cell wall biosynthesis stages ( Figure 3). Interestingly, the direction of expression bias changed between developmental stages in these three genes. For example, expression of homeolog pair was biased in favor of the D T subgenome for Gorai.002G223800 at initiation (1 DPA), but switched to favor the A T subgenome at elongation (5-16 DPA) in wild type developing fibers, whereas expression was biased in favor of the D T subgenome across all evaluated time points in mutant fibers. These results demonstrate that the Li 2 mutation had a greater effect on the D T subgenome and also influenced direction of expression bias for some genes across developmental stages.

Mutation Effects on Functional Distribution of Homeologous Genes during Fiber Elongation
The greater effect of Li 2 on D T biased genes was observed in overall transcript data. In general, the subgenomes contributed unequally to different biological processes [16]; therefore diverse mutation effects could be expected on different functional categories of genes. To determine which biological processes were affected by the mutation, MapMan ontology was used (Data S2).
The distribution of genes from the A T and D T subgenomes with significantly changed expression levels in the mutant were categorized into MapMan functional categories ( Figure 4; Table  S2 in File S1). Relative gene frequencies in functional categories were represented in percents of biased genes in each subgenome (2,958 A T biased genes and 1,620 D T biased genes). Most functional categories were biased in favor of the D T subgenome with the exception of photosynthesis and redox, which only contained A T homeologs. Two functional categories were significantly (p-value ,0.05; Fisher's exact test) biased in enrichment among D T biased genes: secondary metabolism and stress (Table  S2 in File S1). These results demonstrate that different biological processes were unequally affected by Li 2 mutation.
Transcriptional regulators. Transcriptional regulators (TRs) were identified in the G. raimondii genome based on similarity to Arabidopsis TRs and categorized into 76 families. Among them, 229 homeolog pairs were A T biased and 111 were D T biased in elongating cotton fibers. Of the 229 A T biased TRs, 21 (9.2%) of them changed transcription level, whereas of the 111 D T biased TRs 14 (12.6%) of them changed transcription level (Table 3), but this difference was statistically insignificant. Expression levels for the majority of subgenome biased homeologs decreased as the result of Li 2 mutation. However, six TRs (including both homeologs) had increased expression levels in mutant fibers. Three classes of TRs were the most abundant, including Aux/IAA (6 members), bHLH (5 members) and MYB (3 members). Interestingly, two of the three members of MYB TRs had increased expression due to mutation.
Cell wall. In the cell wall functional category, 60 homeologs were A T biased and 40 were D T biased in fiber transcriptome. Ten (16.7%) of the A T biased homeologs changed expression levels; whereas 12 (30%) of the D T biased homeologs changed expression levels (Table 4), indication a higher, but statistically insignificant effect on D T biased homeologs. Interestingly, more D T homeologs (11) than A T homeologs (4) increased transcript levels as a result of the Li 2 mutation. Genes encoding enzymes involved in polysac- charide degradation (14 genes) and cell wall proteins (9 genes) were the most abundant classes.

Validation of Illumina RNA-seq Expression and Subgenome Specific Categorization of Reads by RT-qPCR Analysis
To test the reliability of Illumina sequencing and SNP-based categorization of reads to the A T or D T subgenome of allotetraploid G. hirsutum, RT-qPCR analysis was performed for a subset of 8 genes (selected from Table S1 in File S1) expressed only in WT or Li 2 NILs, and for a subset of 11 genes (selected from Tables 3 and 4) that showed subgenome biased expression. Overall, the results of RT-qPCR analysis were consistent with results of RNA-seq analysis for 19 selected genes (Figures S1 and S2). RT-qPCR analysis confirmed silencing or activation of the expression by the Li 2 mutation for the subset of 8 genes ( Figure  S1). Correlation analysis of the expression patterns revealed strong correlations between RNA-seq and RT-qPCR data. In the subset of 11subgenome biased genes, 7 genes showed 100% correlation (p-value ,0.05) and 4 genes showed 99% correlation (p-value . 0.05; Figure S2).

Discussion
Our results demonstrate that the A T subgenome in general contributed approximately two times more significantly induced genes to the fiber transcriptome than the D T subgenome; however, the Li 2 mutation had greater effects on the D T subgenome than the A T subgenome.

Global Transcript Changes in Subgenomes of G. hirsutum Following Li 2 Mutation
The role of the A T and D T subgenomes in determination of fiber quality in allotetraploid cotton has been extensively discussed in the literature. Allopolyploidization resulted in significant improvements in the desirable agronomic fiber traits in the allotetraploid species in comparison with the diploid progenitors [17,18]. The first evidences showing that QTLs for fiber quality (including length, strength and fineness) were associated with DNA markers mapped to the D T subgenome rather than the A T subgenome was published by Paterson's group [18]. Review of numerous QTLs published from 1998 to 2007 confirmed the observation that the D T subgenome plays a larger role in genetic control of fiber traits [19]. A microarray study published by Wendel's group found that the homeolog expression in G. hirsutum was biased in favor of the D T subgenome in fiber cells [16]. Similar results were reported by Lacape and coauthors utilizing deep sequencing approach to analyze the fiber transcriptome of two allotetraploid species G. hirsutum and G. barbadense [20]. From an evolutionary point of view, these observations are surprising since the genes responsible for improved fiber properties evolved in the diploid AA genome before polyploidization [21]. None of the DD genome species produce spinnable fibers [22].
There are discrepancies in the literature regarding homeolog bias in contribution to fiber traits. Using a core set of 111 RFLP markers, Ulloa and coauthors revealed that the A T subgenome exhibited 68% of QTLs from the five chromosomes, whereas the D T subgenome exhibited only 32% of QTLs from the three chromosomes [23]. Another study utilizing combinations of markers found more fiber trait QTLs in the A T subgenome than in the D T subgenome [24]. The expression analysis of ESTs derived from immature ovules of G. hirsutum TM-1 revealed significant enrichment in all functional categories for A T subgenome ESTs [25].   These inconsistencies could be explained by technical limitations. The QTL studies reported in the current literature are detecting only a small subset of the genes related to fiber traits that may not cover the whole genome and could be insufficient to conclude which subgenome more significantly contributes to fiber properties [19]. The microarray studies evaluated a limited number of homeologous gene pairs, resulting in limited statistical power [16,26]. Lacape and coauthors used next generation DNA sequencing technology for fiber transcriptome analysis; however, they evaluated only 617,000 good quality reads from four libraries without biological replication [20]. Unlike previous studies we obtained ,160 million reads per sample for each of two biological replicates (Table 1), providing ,5.6 times coverage of the G. hirsutum genome (,2.83 Gb per haploid [27]), which is more than enough to deliver statistically powerful transcriptional analysis. Our observation of higher expression of A T than D T genes in the fiber transcriptome is consistent with the results of cotton ovules ESTs analysis [25] and reflects the evolutionary role of the AA diploid progenitor in fiber quality traits of allotetraploid cotton.
It is interesting to note that the Li 2 mutation coincides with an increase in the number of expressed genes, but the homeolog expression bias was significantly decreased in Li 2 fiber. How expression of homeologous genes is regulated in polyploids is still unclear, although it could involve altered regulatory interactions and rapid genetic and epigenetic changes in subgenomes [28]. The evolution homeolog-specific expression after polyploidization has been extensively studied in allotetraploid cotton. Higher rates of homeolog expression bias in natural allotetraploids than in hybrid  Table S2 in File S1 provides Fisher's exact test results. MapMan BIN structure was used for functional categorization of genes regulated by Li 2 mutation. Only functional categories with more than 0.06% gene representation are shown here. Carbohydrates combine 6 BIN classes, including major and minor carbohydrates, glycolysis, fermentation, gluconeogenesis and oxidative pentose phosphate pathway. doi:10.1371/journal.pone.0090830.g004 and synthetic polyploid cottons suggested that the extent of homeolog expression bias increases over time from hybridization through evolution [26,29,30]. The Li 2 mutation is negative for desirable fiber quality traits, resulting in extremely short lint fiber. Significant reduction of homeolog expression bias in short fiber suggests that the extent of homeolog expression bias is also important for fiber quality characteristics.
We observed a reciprocal switch for some genes in expression bias between homeologs during fiber developmental stages in the mutant. A high degree of expression differences between homeologous genes that are developmentally and stress regulated was reported in cotton [17,31,32]. A high-resolution genome-specific study of expression profiling for 63 gene pairs in 24 tissues in allopolyploid and their diploid progenitor cotton species demonstrated that the majority of expression differences between homeologs are caused by cis-regulatory divergence between the diploid progenitors; however, some degree of transcriptional neofunctionalization was detected as well [32].
The Li 2 mutation was mapped to the D T subgenome [6][7][8]; however, the mutated gene and the nature of mutation are currently unknown. The greater mutation effect on the D T than on the A T subgenome observed here suggests two possible mechanisms. The network of regulatory interactions may have been interrupted by a mutation in the D T subgenome resulting in transactivation or repression of individual gene expression levels and expression cascades. Alternatively an epigenetic modulation may preferentially target the D T subgenome. It has been shown that small RNAs can control gene expression and epigenetic regulation in response to hybridization [33][34][35][36]. For example, miRNAs in allopolyploid Arabidopsis triggered unequal degradation of parental target genes [33]. Similarly, in rice hybrids small RNA populations inherited from parents were responsible for biased expression [36]. Additional investigations of epigenetic and chromatin level modifications will provide insights into causes of gene expression variation between subgenomes.

TRs and Cell Wall Functional Categories of Genes Regulated by Li 2 Mutation
Previous transcriptomics and metabolomics studies have shown that the Li 2 mutation terminated the cotton fiber elongation process [6,10]; therefore, genes with changed expression level in the mutant could be involved in elongation. In the present work, we described in detail TRs and cell wall functional categories, which are critical for fiber developmental processes. Many genes in this list (Table 3 and Table 4) were not functionally characterized in cotton; although, based on sequence similarity to genes characterized in Arabidopsis, they could be involved with fiber elongation and represent candidates for further functional analysis in cotton. Transcriptional regulators. Many TRs regulated by Li 2 mutation are involved in hormonal signaling and development. Particularly, two AP2/EREBPs and six Aux/IAAs were in the pool of TRs affected by Li 2 mutation (Table 3). Plant hormones are important for fiber development. It is well documented that exogenous applications of auxins and gibberellic acid stimulate the differentiation of fibers and promote elongation, while abscisic acid and cytokinins inhibit fiber growth in an in vitro cotton ovule culture system [37,38]. Among the auxin responsive genes, Gorai.009G132300 and Gorai.010G227800, whose transcript abundances were significantly reduced in the D T genome of Li 2 , showed sequence similarity to grapevine VvIAA19 regulator [39]. Transgenic Arabidopsis plants over expressing VvIAA19 exhibited faster growth, including root elongation and floral transition, than the control, suggesting that grape Aux/ IAA19 protein is likely to play a crucial role as a plant growth regulator. In the group of bHLH family of TRs, Gorai.007G005700 transcript abundance was significantly reduced in the D T genome of Li 2 and showed sequence similarity to Arabidopsis BEE3, one of several redundant positive regulators of brassinosteroids signaling required for normal growth and development [40].
The actin cytoskeleton plays an important role in cell morphogenesis; down-regulation of GhACT1 disrupted the actin cytoskeleton network in fibers that resulted in inhibition of fiber elongation [41]. A GATA type TR Gorai.005G230900, a homolog of Arabidopsis WLIM1, was down-regulated in the D T subgenome of Li 2 ; a recent study revealed that plant LIM-domain containing proteins (LIMs) define a highly specialized actin binding protein family, which contributes to the regulation of actin bundling in virtually all plant cells [42].
Cell wall. The plant cell wall has a dual role during elongation: to sustain the large mechanical forces caused by cell turgor and to permit controlled polymer extension generating more space for protoplast enlargement [43]. The active biosynthesis of matrix polysaccharides along with increased activity of cell wall loosening enzymes has been considered to be associated with cell wall extension [44][45][46][47][48]. Expression levels of genes encoding enzymes involved in xyloglucan and glucuronoxylan biosynthesis were decreased as a result of Li 2 mutation. Particularly, xyloglucan b-galactosyltransferase (Arabidopsis homolog, MUR3 [49]) and xylosyltransferase (IRX9 [50]) were down-regulated in the A T or D T subgenomes of mutant fibers (Table 4).
Among cell wall proteins arabinogalactans were the most abundant members. Arabinogalactan-proteins have been implicated in many processes involved in plant growth and development, including cell expansion [51,52].
Primary cell wall expansibility and strength is in part mediated by a group of enzymes that comprise a large family of cell wall modifying proteins, the xyloglucan endotransglycosylase/hydrolases (XTHs). XTHs are apoplast-localized enzymes that cleave and reattach xyloglucan polymers [53,54]. The role of XTHs in cotton fiber elongation has been demonstrated: transgenic overexpression of GhXTH1 in cotton increased fiber length up to 20% [55]. D T biased Gorai.007G057400 corresponding to GhXTH1 was down-regulated in mutant fiber.

Conclusion
Repeated polyploidization over evolutionary time has played a significant role in adding genetic variation to the genomes of plant species. The evolution of the homeolog expression after polyploidization has been extensively studied in cotton comparing expression profiling between parental diploids and natural and synthetic allopolyploid species. This is the first report that explored the effects of a single mutation on the homeolog expression of allotetraploid cotton. Our results showed that significant reduction of the homeolog expression bias in mutant fiber correlates with negative fiber traits, indicating that the extent of homeolog expression bias is important for fiber quality characteristics. In addition, we observed significantly greater mutation effects on the D T than on the A T subgenome that might be explained by localization of the mutated gene. Additional studies using numerous naturally occurring cotton fiber mutations are needed to confirm these observations. This work will lead to an understanding of how gene regulation between A T and D T homeologs contributes to enhanced fiber morphology in cultivated cotton allopolyploids.

Plant Material and RNA Isolation
The cotton short fiber mutant Li 2 was developed as a nearisogenic line (NIL) with the WT upland cotton line DP5690 as described before [6]. Growth conditions and fiber sampling were previously described [6]. Cotton bolls were harvested at the following time-points during development: day of anthesis (DOA), 1, 3, 5, 8, 12, 16, and 20 days post-anthesis (DPA). Cotton fibers were isolated from developing ovules using a glass bead shearing technique to separate fibers from the ovules [56]. Total RNA was isolated from detached fibers using the Sigma Spectrum Plant Total RNA Kit (Sigma-Aldrich, St. Louis, MO) with the optional on column DNase1 digestion according to the manufacturer's protocol. The concentration of each RNA sample was determined using a NanoDrop 2000 spectrophotometer (NanoDrop Technologies Inc., Wilmington, DE). The RNA quality for each sample was determined by RNA integrity number (RIN) using an Agilent Bioanalyzer 2100 and the RNA 6000 Nano Kit Chip (Agilent Technologies Inc., Santa Clara, CA) with 250 ng of total RNA per sample.

RT-qPCR Analysis
The experimental procedures and data analysis related to RT-qPCR were performed according to the Minimum Information for Publication of Quantitative Real-Time PCR Experiments (MIQE) guidelines [57]. Eight fiber developmental time-points mentioned above were used for RT-qPCR analyses of homeolog pairs which showed reciprocal expression biases. Only one time point, 8 DPA, was used for RT-qPCR confirmation of RNA-seq data of selected genes. The detailed description of reverse transcription, qPCR and calculation were previously reported [6]. Single nucleotide polymorphisms that distinguish the A T and D T subgenome copies of the selected genes were identified by aligning reads from the RNA-seq data to the G. raimondii reference mRNA sequences [11]. These homeologous SNPs were used to design subgenome specific primers by the SNAPER approach, whereby an additional mismatch is included near the end of the SNP-specific primers to increase stringency [58]. Primer sequences are provided in Table S3 and Table S4 in File S1. Correlations of biased expression patterns between RNA-seq and RT-qPCR data were calculated using GraphPad Prism 5 software (Pearson test).

Library Preparation and Sequencing
RNA samples from Li 2 and wild type cotton fiber at 8 DPA (in two biological replicates) were subjected to paired-end Illumina mRNA sequencing (RNA-seq). Library preparation and sequencing were conducted by Data2Bio LLC (2079 Roy J. Carver Co-Laboratory, Ames, Iowa). Indexed libraries were prepared using the Illumina protocol outlined in the TruSeq RNA Sample Prep Guide (Part# 15008136 Rev. A, November 2010). The library size and concentration were determined using an Agilent Bioanalyzer. The indexed libraries were combined and seeded onto one lane of the flowcell. The libraries were sequenced using 101cycles of chemistry and imaging, resulting in paired end (PE) sequencing reads with length of 26101 bp. The raw reads were submitted to the Sequence Read Archive (accession number SRP026301).

Processing of Illumina RNA-Seq Reads and Mapping to A T and D T Subgenomes of Gossypium hirsutum
The reads were trimmed with SICKLE (https://github.com/ najoshi/sickle) using a quality score cutoff of 20. Mapping the reads (in pairs where both reads of a pair passed trimming) to the 13 chromosomes of the G. raimondii genome D 5 v2 reference sequence was performed using GSNAP [59]. Default parameters were used, but with the flags ''-n 1 -Q'' which means that only a single mapping was reported for each read, and reads with multiple equally good hits were thrown away rather than randomly mapped. We used a cotton SNP index generated between DD genome G. raimondii and the AA genome G. arboreum to categorize reads of the allotetraploid G. hirsutum as belonging to the A T or D T subgenomes according to the method reported previously [14].

Digital Gene Expression Analysis
The comparison of the number of reads mapped to the genes of G. raimondii reference genome was used as an indicator of the relative digital gene expression (DGE). The JMP/Genomics 6.0 (SAS, Cary, NC, USA) was used for data normalization and statistical analysis. The data was normalized using TMM (Trimmed Mean of M component) method [60]. Genes with less than 10 reads in one sample were removed before normalization; from 37,223 genes assigned to chromosomes, 31,114 genes passed filtering conditions and were processed for normalization. The ANOVA process was fit to the normalized data, with the data following a Poisson distribution. This was accomplished with a generalized linear mixed model for each gene: Y ij = T i + G j + TG ij + E ijk , where T is the treatment effect for the ith biological treatment (Li 2 or wild-type fiber), G is the specific subgenome type effect for the jth subgenome type (A T , D T , X and N categorized reads), their interaction (TG), and the error term (E). The linear model was used to test the null hypothesis that expression of a given gene was not different. Specifically, multiple comparisons were made between fiber type (Li 2 versus wild type) and A T /D T subgenomes as well as combinations of these factors, such as fiber type in A T and D T subgenomes. We identified genes for which the difference in expression levels within these a priori questions were significantly different (false discovery rate#0.05) [61].

Functional Categorization of Genes
Functional categorization of genes was performed using Map-Man ontology [62]; the MapMan mapping for G. raimondii is available at http://mapman.gabipd.org/. Fisher's exact test was used to estimate enrichment or depletion relative to background of functional categories with differentially regulated genes. Figure S1 RT-qPCR confirmation of silencing or activation of genes as a result of mutation. Bar charts represent RNA-seq and RT-qPCR data (side by side) at 8 DPA of fiber development for 8 randomly selected genes from Table S1 in File S1. Error bars indicate standard deviation from two biological replicates for RNA-seq data and three biological replicates for RT-qPCR. (TIF) Figure S2 RT-qPCR confirmation of biased expression of homeolog pairs. Bar charts represent RNA-seq and RT-qPCR data (side by side) at 8 DPA of fiber development for 11 randomly selected genes from Table 3 and Table 4. Pearson correlation (GraphPad Prism 5 software) of expression patterns for selected genes between RNA-seq and RT-qPCR data is provided in the table; correlation coefficients with p-value less than 0.05 are shown in boldface and underlined. Error bars indicate standard deviation from two biological replicates for RNA-seq data and three biological replicates for RT-qPCR. (TIF)

Supporting Information
File S1 Supporting tables. Table S1. Silencing or activation of genes as a result of mutation. Table S2. Mutation effects on functional distribution of homeolog genes. Fisher's exact test results. Table S3. Primer's sequences for detection expression of homeolog pairs. Table S4. Primer's sequences. (DOCX) Data S1 Statistical data for significantly regulated genes.

(TXT)
Data S2 A T /D T biased genes annotated by MapMan ontology. (TXT)


extfieldsofstudy:
Materials Science
Medicine

text:
Mapping Photothermally Induced Gene Expression in Living Cells and Tissues by Nanorod-Locked Nucleic Acid Complexes

The photothermal effect of plasmonic nanostructures has numerous applications, such as cancer therapy, photonic gene circuit, large cargo delivery, and nanostructure-enhanced laser tweezers. The photothermal operation can also induce unwanted physical and biochemical effects, which potentially alter the cell behaviors. However, there is a lack of techniques for characterizing the dynamic cell responses near the site of photothermal operation with high spatiotemporal resolution. In this work, we show that the incorporation of locked nucleic acid probes with gold nanorods allows photothermal manipulation and real-time monitoring of gene expression near the area of irradiation in living cells and animal tissues. The multimodal gold nanorod serves as an endocytic delivery reagent to transport the probes into the cells, a fluorescence quencher and a binding competitor to detect intracellular mRNA, and a plasmonic photothermal transducer to induce cell ablation. We demonstrate the ability of the gold nanorod-locked nucleic acid complex for detecting the spatiotemporal gene expression in viable cells and tissues and inducing photothermal ablation of single cells. Using the gold nanorod-locked nucleic acid complex, we systematically characterize the dynamic cellular heat shock responses near the site of photothermal operation. The gold nanorod-locked nucleic acid complex enables mapping of intracellular gene expressions and analyzes the photothermal effects of nanostructures toward various biomedical applications.

Dividing the number of GNRs determined by UV-Vis extinction measurement yield the number of LNA probe attached to a single GNR to approximately 34 LNA per GNR.
Considering each GNR carries ~34 copies of LNA probes in a cell, the final probe concentration inside a cell was approximately 45 nM. For cell imaging, 2×10 11 GNR/ml GNR-LNA complexes targeting heat shock protein 70 (HPS70) mRNA, β-actin mRNA, and random sequence were tested.

Quantification of internalized GNRs
Inductively coupled plasma mass spectrometry (ICP-MS) was applied to quantify GNRs internalizations. Briefly cells were cultured overnight. GNRs with concentration of 10 10 , 2×10 11 , and 5×10 11 GNR/ml were introduced to 24 well plates and incubated for 12 hours. Cells were then washed for 3 times with 1X PBS to remove excessive GNRs that were not internalized. Cells were harvested with 0.25% trypsin solution and centrifuged at 1500 rpm for 5 min. The cell pellets were then resuspended in 1X PBS to obtain a final concentration of 5×10 5 cell/ml. The GNR samples were prepared by acid digestion of cells using the CEM MARS6 microwave digestion system (Matthews, N.C.).
The sample was weighed in a tared Teflon vessel. Freshly-prepared aqua regia (Aristar-Plus grade HNO3:HCl at 1:3 ratio) was added. After one hour of pre-digestion, 1.0 ml of 18 megohm water was added prior to sealing the vessels. During the digestion process, the temperature of the solution reached to 200 o C after 20 min, and was held at this temperature for an additional 15 min. After cooling, the digest was transferred into a pre-weighed 15mL trace-metal-free PP tube (VWR Scientific). The mass of GNRs in the solution were then determined by Agilent 7700x CIP-MS system (Agilent Technologies, Santa Clara, CA) and converted to the number of GNRs per cell using the calculated number of gold atom in each nanorods according to the equation previously reported 10 .

Statistical analysis.
All experiments were performed independently for at least three times. For mapping the spatial and temporal heat shock response of cells after laser irradiation, at least 50 cells were measured for each data point. The data are expressed as mean ± SEM. For characterizing laser induced temperature elevation, the experiments were repeated (n=3) for statistical analysis. Statistical analysis was carried out using Student's t-test.

Temperature measurement
The temperature was measured at different time points and locations from the laser.
For thermal heating and cooling experiments, miniaturized thermocouples (Omega) with 10 µm diameter were placed inside the media near the cell monolayer. Data were recorded using a thermocouple reader (Okolab) at different time points. The thermocouple was placed at different locations using a costumed 3D traversing system to determine the temperature distribution near the laser. To verify the temperature measurement and resolve the spatial temperature distribution, infrared thermometry was also performed using a mid-wavelength infrared camera (SC6700, FLIR Thermovision

B. Numerical analysis of temperature profile Heat equation
To investigate the temperature distribution near the site of photothermal operation and the effects of laser power and duration, we performed computational analyses to estimate the temperature profile near the irradiation area. In particular, the heat diffusion equation, which describes the distribution of temperature over time, was considered in the analysis 11 . A heating term was introduced to represent the localized heat generated due to the photothermal effect (equation S1).
where T is the temperature distribution and S is the local photothermal heating source.
D is the thermal diffusivity of the medium. Thermal diffusivity of water was used in the calculation. The value can be modified to describe the heat transfer in other media.
The localized heating is modeled by assigning an appropriate temperature at the location of photothermal operation. The partial differential equation was written in twodimensional form to approximate the photothermal ablation experiment (equation S2).
The equation was discretized in space and time using the finite difference approach and implemented in MATLAB (equation S3). 12 n, i, and j are discrete time and space indexes.

Boundary and initial conditions
Supplementary Fig. S4a shows the domain considered in the simulation. In the simulation domain, there were 20 by 80 elements and each element was 100 µm by 100 µm in size. Initial and boundary conditions were required to solve the heat equation.
For initial condition, the temperature was considered to be uniform at the reference temperature initially. The thermal boundary condition was determined by the heat transfer through the cover slip ( Supplementary Fig. S4b). In particular, the heat loss to the cover glass was determined by the heat transfer coefficients at the liquid-glass interface and glass-air interface, and conductive heat transfer through the cover glass.
T f , T s , T b , and T ∞ were the temperature of the liquid, the liquid-glass interface, the glassair interface, and the ambient. In our experimental condition, the heat flux was mainly determined by the heat transfer coefficient at the glass-air interface. Since most of the heat loss occurred near the laser source on the cover glass, the heat transfer was considered to be the same in all boundaries. The parameters and values used in the numerical simulation are described in Table S2. 11

Temperature distribution with continuous laser irradiation
The temperature profiles with continuous laser heating were shown in Fig. 2h in the main text (see also Supplementary movie S3). In the simulation, heat diffused rapidly to the surrounding media in the first few seconds. The temperature profile then gradually reached a quasi-steady state value when the heat flux is approximately equal to the heat loss to the substrate. Supplementary Fig. S5 shows the temperature profiles near the laser. The region of temperature elevation was localized near the heat source. The temperature dropped to 50% of the maximum value in 200-300 µm, representing a characteristic heat transfer length scale of the system. This length scale was in good agreement with the value measured in the experiment by infrared thermometry and miniaturized thermocouples (see Fig. 2e).

Heating and cooling time constants
To estimate the thermal time constants

Effects of laser duration
We then investigated the effects of the duration of heating on the temperature profile.

C.
Supplementary Tables   Table S1. Sequences of single-stranded LNA probes with alternating DNA/LNA monomers for intracellular detection. Bold italic letters represent LNA monomers.