text:
Neuroprotective effect of nerolidol against neuroinflammation and oxidative stress induced by rotenone

Background Parkinson disease (PD) is a movement disorder affecting 1 % of people over the age of 60. The etiology of the disease is unknown; however, accumulating evidence suggests that mitochondrial defects, oxidative stress, and neuroinflammation play important roles in developing the disease. Current medications for PD can only improve its symptoms, but are unable to halt its progressive nature. Although many therapeutic approaches are available, new drugs are urgently needed for the treatment of PD. Thus, the present study was undertaken to investigate the neuroprotective potential of nerolidol, a sesquiterpene alcohol, on a rotenone-induced experimental model of PD, where male Wistar rats intraperitoneally received rotenone (ROT) at a dose of 2.5 mg/kg of body weight once daily for 4 weeks. Results Nerolidol, which has antioxidant and anti-inflammatory properties, was injected intraperitoneally at 50 mg/kg of body weight, once daily for 4 weeks, and at 30 min prior to ROT administration. ROT administration significantly reduced the activities of antioxidant enzymes superoxide dismutase (SOD) and catalase (CAT), and the level of the antioxidant tripeptide glutathione (GSH). Moreover, ROT increased the levels of the lipid peroxidation product malondialdehyde (MDA), proinflammatory cytokines (IL-1β, IL-6, and TNF-α), and inflammatory mediators (COX-2 and iNOS) in rat brain tissues. Immunostaining of brain tissue sections revealed a significant increase in the number of activated astrocytes (GFAP) and microglia (Iba-1), along with the concomitant loss of dopamine (DA) neurons in the substantia nigra pars compacta and dopaminergic nerve fibers in the striatum of ROT-treated rats. As expected, nerolidol supplementation to ROT-injected rats significantly increased the level of SOD, CAT, and GSH, and decreased the level of MDA. Nerolidol also inhibited the release of proinflammatory cytokines and inflammatory mediators. Finally, nerolidol treatment prevented ROT-induced glial cell activation and the loss of dopaminergic neurons and nerve fibers, and ultimately attenuated ROT-induced dopaminergic neurodegeneration. Conclusion Our findings are the first to show that the neuroprotective effect of nerolidol is mediated through its anti-oxidant and anti-inflammatory activities, which strongly supports its therapeutic potential for the treatment of PD.

Background
Parkinson disease is a neurodegenerative disorder that is characterized by a progressive loss of dopamine (DA) neurons in the substantia nigra pars compacta (SNc) and dopaminergic nerve terminal fibers of the striatum that are enriched in the DA transporter [1]. As a result of dopaminergic neurodegeneration, the level of the DA neurotransmitter (which mediates the control of movement) is significantly reduced in the striatum. Therefore, dopaminergic neurodegeneration leads to impairments in movement, a resting tremor, rigidity, and a disturbance in gait. Although the etiology for dopaminergic neuronal degeneration in PD is not completely understood, current evidence suggests that oxidative stress and neuroinflammation play key roles in the pathogenesis of PD [2,3].
Rotenone (ROT), an inhibitor of mitochondrial complex I, has been widely used as a herbicide/pesticide in agriculture. A ROT challenge, due to its complex I inhibitory property, induces pathological features in animals that are similar to those seen in human PD patients. Thus, ROT-treated animals represent a promising animal model with construct validity [4][5][6]. The ROT model recapitulates most of the pathological features observed in human PD pathogenesis, including the loss of DA neurons in SNc, and enhanced oxidative stress and neuroinflammation in the nigrostriatal dopaminergic pathway [2,[6][7][8]. Moreover, a ROT challenge to rodents induced the formation of α-synuclein cytoplasmic inclusions in DA neurons, Lewy pathology, DJ-1 acidification and translocation, proteasomal dysfunction, and nigral iron accumulation [9]. Therefore, the ROT-treated rat is a suitable model for investigating novel therapeutic agents targeting oxidative stress and neuroinflammation in PD.
Currently available drugs alleviate the symptoms of PD, but are inadequate in terms of halting the progression of the disease. Additionally, motor complications during advanced stages of the disease, the adverse effects of the available drugs, and non-motor symptoms remain huge challenges for long-term therapy. Therefore, newer therapeutic agents and approaches are urgently needed to stop or delay the progressive nature of the disease [10]. PD is intimately connected to excessive oxidative stress that damages brain areas and accelerates the process of neurodegeneration, especially in the nigrostriatal area. This process is an important breakpoint for controlling the disease. Consequently, in addition to other therapeutic approaches, treatment with antioxidants is gradually gaining importance in the pharmacotherapy of PD.
Plants are a potential source of novel drugs to benefit mankind. Much research effort has focused on discovering new antioxidant compounds from plants. An increasing number of studies have shown that essential oils obtained from medicinal plants exhibit a variety of biological properties, such as anticonvulsant [11], analgesic [12], anxiolytic [13], antidepressant [14], antioxidant [15], and anti-inflammatory [16] properties. Some of these effects are frequently attributed to terpenes, which are chemical components that are present in these essential oils [11,17]. Nerolidol (NRD) [3,7,11-trimethyl-1,6,10dodecatrien-3-ol] is a sesquiterpene alcohol (Fig. 1) that is found in essential oils from plants such as Baccharis dracunculifolia [18], Amaranthus retroflexus [19], and Canarium schweinfurthii [20]. Nerolidol is an aliphatic sesquiterpene that is commonly used in cosmetics, perfumes, shampoos, and soaps [21], as well as non-cosmetic products (cleansers and detergents). Essential oils from Ferula fukanensis containing nerolidol have been shown to prevent nitric oxide (NO) production and NO-induced gene expression, suggesting its plausible use as an antioxidant agent. NRD is hydrophobic in nature with an XlogP3 value of 4.6 (>2.5 is needed for efficient transport across blood brain barrier), and thus NRD can cross easily the blood brain barrier. NRD significantly reduces the level of lipid peroxidation and nitrite content in the hippocampus of mice, which protects against oxidative stress [15], suggesting an important antioxidant property for NRD. A recent study reported that NRD reduces interleukin-1β (IL-1β) production in lipopolysaccharide-induced peritoneal macrophages [16]. Additionally, NRD shows potent anti-inflammatory activity through the suppression of IL-1β and tumor necrosis factor-alpha (TNF-α) levels in an experimental mouse model of pain [16]. Therefore, the main objective of the present study was to investigate the anti-oxidative and anti-inflammatory effects of NRD on ROT-induced neurodegeneration in rats.

Drugs and chemicals
Polyclonal rabbit anti-cyclooxygenase-2 (COX-2), antiinducible nitric oxide synthase (iNOS), and anti-glial fibrillary acidic protein (GFAP) antibodies were purchased from Abcam, Cambridge, MA, USA. Anti-ionized calcium-binding adaptor molecule-1 (Iba-1) polyclonal rabbit antibody was procured from Wako Chemicals, Richmond, VA, USA. Polyclonal rabbit anti-tyrosine hydroxylase antibody was obtained from Novus Biologicals, Littleton, CO, USA. Alexa fluor 488 conjugated secondary goat anti-rabbit antibodies were purchased from Life Technologies, Grand Island, NY, USA. Biotinylated secondary anti-rabbit antibody was purchased from Jackson Immunoresearch, West Grove, PA, USA. Nerolidol was purchased from Santa Cruz Biotechnology Inc., CA, USA. ROT and the assay kit for reduced glutathione (GSH) and other reagents of analytical grade were purchased from Sigma-Aldrich, St. Louis, MO, USA.

Experimental animals
We used male Wistar rats at 6-7 months old (weighing 280-300 g) from the Animal Research Facility of the College of Medicine and Health Sciences, United Arab Emirates University, UAE. Prior to the start of the experiment, a maximum of 4 rats were housed per cage, and rats were acclimatized for 1 week to the laboratory conditions. The cages of the animals were changed twice a week. The animals were housed and kept under standard laboratory light and dark cycle conditions. The animals were fed with commercially available rodent food and water ad libitum. All experiments were conducted between 0900 and 1500 hours. The Animal Ethics Committee of the United Arab Emirates University, UAE, approved the experimental protocol for animal experimentation.

Experimental design
To induce the PD in rats, ROT (2.5 mg/kg body weight) was administered intraperitoneally (i.p.) once daily for 4 weeks. The regimen used in the current study for the induction of Parkinsonism in rats was adopted as reported earlier [22]. Briefly, ROT was first dissolved in dimethyl sulfoxide (DMSO) to prepare a 50X stock solution and stored at −80 °C for further use. Before injection, the stock solution of ROT was thawed and further diluted in sunflower oil to obtain a final concentration of 2.5 mg/mL. To evaluate its neuroprotective efficacy, NRD was initially diluted in olive oil to obtain a final concentration of 50 mg/2 mL, and injected i.p. 30 min prior to ROT administration at a dose of 50 mg/kg of body weight once daily for 4 weeks. The dose of NRD was selected based on dose-response studies (data not shown) as well as earlier reports [15,23]. Rats that were used as controls received an equivalent volume of vehicle only. The animals were sacrificed 48 h after the last injection of NRD or ROT or both in combination to completely eliminate these drugs from the body. The rats were divided into 4 experimental groups where each group consisted of 8 animals; the groups were named as follows: Group I, the vehicle-injected control group (CONT). Group II, the rotenone-injected and vehicle-treated group (ROT). Group III, the rotenone-injected and NRD-treated group (ROT + NRD). Group IV, the NRD-only injected group (NRD).

Tissue collection
At the end of the experiment, the animals from all groups were anaesthetized with pentobarbital (40 mg/kg of body weight) and a cardiac perfusion was performed using 0.01 M of phosphate-buffered saline (PBS) pH 7.4 to completely clear the blood. The skull of the rat was open to quickly isolate the brain. The brain was placed on an ice-plate and cut along the midline to separate the 2 cerebral hemispheres. For biochemical assay, the midbrain and striatum regions of the brain were isolated from 1 hemisphere and immediately fresh frozen in liquid nitrogen for further use. The other hemisphere of the brain was incubated in 4 % paraformaldehyde solution for 48 h and subsequently exchanged with 10 % sucrose solution containing 0.1 M of phosphate buffer (PB) 3 times daily for 3 consecutive days at 4 °C prior to cryostat sectioning.

Sample preparation for biochemical assay
For biochemical assay, the midbrain tissue of each animal was homogenized in KCl buffer (10-mM Tris-HCl, 140-mM-NaCl, 300-mM KCl, 1-mM EDTA, and 0.5 % Triton X-100) at pH 8.0 supplemented with protease and phosphatase inhibitor, keeping the samples on ice. The tissue homogenates of each sample were centrifuged at 14,000×g for 20 min at 4 °C to obtain the post-mitochondrial supernatant (PMS) fraction. This PMS fraction was used to estimate the levels of the antioxidant enzyme GSH, lipid peroxidation product, and proinflammatory cytokines using spectrophotometric measurements and an enzyme-linked immunosorbent assay (ELISA) following a standard protocol, as reported earlier [22].

Assay for lipid peroxidation
A malondialdehyde (MDA) assay kit procured from Northwest Life Science (Vancouver, WA, USA) was used to determine the amount of lipid peroxidation, as reported earlier [22]. Briefly, the samples or standards (250 µL) were incubated in the presence of acid reagent (250 µL) and thiobarbituric acid (250 µL) and vortexed vigorously. Samples were further incubated for 60 min at 60 °C followed by centrifugation at 10,000×g for 2-3 min. The reaction mixture was then aseptically transferred without disturbing the pellet to a cuvette and the absorbance was recorded at 532 nm. The concentration of MDA was calculated using a standard curve and expressed as µM MDA/mg protein.

Estimation of GSH
A GSH kit was used to estimate the GSH level, as reported earlier [22]. Briefly, the samples were first deproteinized with 5 % 5-sulfosalicylic acid solution and centrifuged to remove the precipitated protein. The supernatant was used to measure the GSH level. Samples and standards of different concentrations (10 µL) were added in each well of a 96-well plate and incubated for 5 min with the working mixture (150 µL; assay buffer + 5,5′-dithiobis (2-nitrobenzoic acid) + glutathione reductase). Diluted NADPH solution (50 µL) was added to each well and mixed properly. The absorbance of the samples was measured at 412 nm with kinetics capability for 5 min using a microplate reader. The results were expressed as µM GSH/mg protein.

Assay for antioxidant enzymes activity
The activity of antioxidant enzymes such as superoxide dismutase (SOD) and catalase (CAT) were determined following the manufacturer's instructions of a kit (Cayman Chemicals Company, Ann Arbor, MI, USA), as reported earlier [22]. The CAT activity was expressed as nmol/min/mg protein and the SOD activity was expressed as U/mg protein.

Assay for pro-inflammatory cytokines
Commercially available ELISA kits for IL-1β, IL-6, and TNF-α were purchased from R&D systems, Minneapolis, MN, USA. The level of IL-1β, IL-6, and TNF-α were estimated as described earlier [22]. The results were expressed as pg/mg protein.

Immunohistochemistry for tyrosine hydroxylase (TH) expression
The one hemisphere of the brains collected from each animal were serially cut after fixation and processed for immunohistochemical analysis. Briefly, 14 μm-thick coronal brain sections were cut at the level of the striatum and SNc using a cryostat (Leica, Wetzlar, Germany) for the immunohistochemical analysis of TH. Sections were washed twice with 0.01 M of PBS, pH 7.4, and then incubated with blocking reagent (10 % normal goat serum in PBS containing 0.3 % Triton-X 100) for 1 h. Next, the sections were incubated overnight with a primary polyclonal rabbit antibody against TH (1:500) at 4 °C. Sections were washed and incubated with biotinylated secondary antirabbit (1:1000) antibody for 1 h at room temperature. The brain sections were incubated with avidin-biotin complex (Vector Laboratories Ltd. Burlingame, CA, USA) and 3,3′ diaminobenzidine (DAB) to visualize and analyze the TH immunoreactivity. Finally, the sections were coverslipped using DPX mounting medium. The slides were then viewed under a light microscope (Olympus, Hamburg, Germany) and images were acquired for analysis.

Immunofluorescence staining of GFAP and Iba-1
Immunofluorescence staining was performed using 14 µm-thick striatum sections to quantify the number of activated GFAP-positive astrocytes and Iba-1-positive microglia. The brain sections were first washed twice with PBS and incubated with blocking reagent (10 % normal goat serum in PBS containing 0.3 % Triton-X 100) for 1 h. The sections were then incubated overnight with the primary polyclonal rabbit antibodies ant-GFAP (1:1000) and anti-Iba-1 (1:1000) at 4 °C. The sections were washed and incubated with fluorescent secondary anti-rabbit Alexa fluor 488 antibody for 1 h at room temperature. Sections were then washed, mounted on slides, and coverslipped using the mounting medium Fluoroshield (Sigma Aldrich, St. Louis, MO, USA). The images were captured using a fluorescence microscope EVOS FL (Thermo Fisher Scientific, Waltham, MA, USA).

Assessment of TH-immunopositive (TH+) dopamine (DA) neurons in the SNc and TH-immunoreactive (TH-ir) DA nerve fibers in the striatum
To evaluate the ROT-induced neurodegeneration and neuroprotective effect of NRD, the total number of TH+ DA neurons at 3 different anatomical levels of the SNc (−4.8, −5.04, and −5.28 mm of the bregma) were counted. The average number of TH+ neurons were calculated and converted as a percentage with reference to the control. The loss of striatal fibers was evaluated by measuring the intensity of TH-ir dopaminergic fibers in the striatum using Image J software (NIH, Bethesda, MD, USA). The intensity of the TH-immunoreactive nerve fibers in 3 different fields of brain sections (3 sections per animal) within the striatal region (adjacent to 0.3 mm of the bregma) was measured to examine the DA nerve fibers loss. An average of the 3 sections was calculated and presented as a percentage with reference to the values of the control group. The intensity of the overlying cortex area was used as the background measurement and subtracted from the value generated from the striatum. An investigator that was 'blind' to the experimental groups counted the TH+ DA neurons and determined the immunoreaction intensity of the TH fibers.

Assessment of activated astrocytes and microglia in the striatum
Three coronal sections from a comparable anatomical level of striatum from each animal were used to analyze and count the number of activated astrocytes and microglia. Activated astrocytes and microglia were considered for counting based on the intense immunoreactivity of GFAP-and Iba-1-labeled cells, whose characteristic morphological features include hypertrophied and extended glial processes. The total number of activated astrocytes and microglia were counted from 3 randomly chosen fields of an equal area in each section using the Image J software (NIH, Bethesda, MD, USA), and the results were presented as a percentage.

Western blot analysis of COX-2 and iNOS expression
Western blot analysis was carried out to measure the level of COX-2 and iNOS expression in different groups of animals following the protocol, as reported earlier [22]. Briefly, striatal tissues isolated from each animal were homogenized in radioimmuno-precipitation buffer with protease and phosphatase inhibitors. The cell lysates were then centrifuged at 15,000 rpm for 20 min. The supernatant containing cytoplasmic fractions was isolated, and protein concentration was measured as mentioned below. The cytoplasmic fraction containing equal amounts of protein (35 μg) were loaded and separated using 10 % SDS-polyacrylamide gel electrophoresis. The proteins were then transferred onto a PVDF membrane and incubated overnight at 4 °C with specific primary rabbit polyclonal antibodies against COX-2 (1:1000) and iNOS (1:500). The membrane was washed and then incubated with horseradish peroxidase-conjugated secondary anti-rabbit antibody. The protein recognized by the antibody was visualized using an enhanced chemiluminescence Pico kit (Thermo Fisher Scientific, Rockford, IL, USA). The blots were stripped and re-probed for β-actin (1:5000; monoclonal mouse, Millipore, MA, USA) as a loading control. The intensity of the bands was measured using densitometry and quantified using Image J software (NIH, Bethesda, USA).

Protein estimation
The concentration of protein in each sample was estimated using the Pierce BCA protein assay kit (Thermo Fisher Scientific, Rockford, IL, USA) following the manufacturer's instructions.

Statistical analyses
The data were expressed as the mean value ± SEM. The data were analyzed with Graph Pad (InStat software, La Jolla, CA, USA) using a one-way analysis of variance (ANOVA) followed by a Tukey's test to determine the statistical significance between various groups. For all of the tests, the criterion for any statistically significant difference was set at p < 0.05.

NRD administration prevents the loss of DA neurons in the SNc and DA nerve fibers in the striatum
To investigate the beneficial effect of NRD on DA neurodegeneration in ROT-injected rats, TH immunohistochemical analysis was performed to assess the expression of healthy TH+ DA neurons in the SNc and the TH-ir DA nerve fibers density in the striatum. The results show that chronic ROT treatment caused a significant loss (p < 0.05) of DA neurons in the SNc area as expected when compared to vehicle-injected control rats (Fig. 2a,  c). In contrast, NRD supplementation prior to ROT injection in rats provided significant (p < 0.05) protection to DA neurons when compared to only ROT injection of rats. The DA neurons of the SNc project their nerve terminals to the striatum where they are highly enriched with the dopamine transporter (DAT). Degeneration of DA neurons in the SNc decreases DAT expression in the striatum, which further confirms the neuronal loss. Therefore, we examined whether the loss of DA neurons in the SNc area is correlated with DA-terminal loss by measuring the intensity of TH-ir dopaminergic nerve terminal fibers of the striatum. A significant decrease (p < 0.05) was observed in the intensity of TH-ir fibers in ROT-injected animals when compared to vehicle-treated controls. However, NRD supplementation prior to ROT administration significantly increased (p < 0.05) TH-ir intensity in nerve terminal fibers, suggesting that NRD has a potent neuroprotective effect against ROT-induced dopaminergic neurodegeneration (Fig. 2b, d).

NRD treatment suppresses excessive lipid peroxidation and ameliorates GSH levels
Next, we tested whether NRD mediated the neuroprotective effect that we observed in our model is due to its antioxidant activity because NRD is an antioxidant. To clarify this issue, the role of NRD in ROT-induced neurodegeneration was investigated by examining the level of MDA, a marker for lipid peroxidation in the midbrain region. As shown in Fig. 3a, ROT treatment resulted in a massive amount of MDA production in the ROT group of animals when compared to control rats. Treatment with NRD significantly abolished (p < 0.05) the elevated MDA production. Increased MDA production due to oxidative stress correlated with a decrease in the bioavailability of the tripeptide antioxidant GSH in midbrain cells. The ROT-challenged rats also showed a significantly decreased (p < 0.05) level of GSH compared to vehicle-injected control rats (Fig. 3b). Supplementation with NRD prior to ROT administration significantly attenuated (p < 0.05) the level of GSH compared to the ROT-injected rats. Taken together, these results suggest that the neuroprotective effect of NRD on ROT-induced neurodegeneration might be mediated through its antioxidative effect.

NRD reverses the ROT-induced decrease in the activity of antioxidant enzymes
We also measured the activity of antioxidant enzymes SOD and CAT. ROT administration significantly decreased (p < 0.01) SOD and CAT activity when compared with the vehicle-injected control group. However, NRD supplementation significantly increased (p < 0.05) the activities of SOD (Fig. 3c) and CAT (Fig. 3d) when compared with the ROT-injected animals. Control animals did not show any significant changes in the activities of SOD and CAT when compared to the NRD-onlyinjected rats.

Neuroprotective effects of NRD involves the inhibition of glial cell activation
Glial cell (astrocyte and microglia) activation following ROT administration has been observed, and is considered to be an index of the inflammatory response. We observed an increase in the expression of GFAP and Iba-1, which are markers of activated astrocytes and microglia, respectively, upon ROT administration.
Immunofluorescence staining showed a significant increase (p < 0.01) in the number of activated astrocytes and microglia in ROT-injected animals compared to the vehicle-injected control animals (Fig. 4a-d). Interestingly, NRD treatment to the ROT-injected rats significantly decreased (p < 0.05) the number of activated astrocytes and microglia compared to the ROT-injected rats. Animals injected with NRD alone did not show any remarkable activation of astrocytes and microglia when compared to control animals.

NRD attenuates the ROT-induced release of proinflammatory cytokines
The persistent activation of astrocytes and microglia accompanied by the sustained secretion of inflammatory mediators is thought to be involved in the neuronal Values are expressed as percent mean ± SEM (n = 3). d A significant decrease (*p < 0.05) in the TH-ir fibers was observed in the ROT group compared to the CONT group. NRD treatment significantly inhibited ( # p < 0.05) the loss of TH-ir fibers in the ROT + NRD group compared to the ROT group. The CONT rats and NRD-only-injected rats did not show a remarkable loss of TH-ir fibers injury in PD. Therefore, the release of proinflammatory cytokines such as IL-1β, IL-6, and TNF-α in response to ROT exposure was examined and quantified using ELISA. A significant increase (p < 0.01) in the IL-1β, IL-6, and TNF-α level was observed in the ROT-injected animals when compared to the vehicle-injected control animals (Fig. 5a-c). However, NRD supplementation to the ROT-administered animals significantly decreased (p < 0.05) the level of these proinflammatory cytokines when compared to the ROT-injected animals (Fig. 5a-c). Animals injected with NRD alone did not show any significant changes in the level of analyzed proinflammatory cytokines when compared to control animals.

Effect of NRD on the expression of COX-2 and iNOS
The expressions of COX-2 and iNOS were also examined using Western blot in tissue lysates that were isolated from the striatum region (Fig. 6a, b). A significant increase (p < 0.05) in COX-2 expression (Fig. 6c) was observed in ROT-treated animals (217.42 vs. 100 %) when compared to the vehicle-injected control group of animals. However, following treatment with NRD in the ROT-administered rats, a remarkable reduction in the level of COX-2 (148.16 vs. 217.42 %) was observed when compared to the ROT-treated rats. ROT administered animals also showed a significant increase (p < 0.05, 206.35 vs. 100 %) in iNOS expression (Fig. 6d)

Discussion
PD is one of the most common progressive neurodegenerative disorders with a complex pathogenesis. The current hypothesis suggests that mitochondrial dysfunction, oxidative stress, and neuroinflammation play a crucial role in the development of PD. Systemic administration of ROT to rats has been shown to induce behavioral and pathological symptoms of PD, such as motor dysfunction, Fig. 3 Quantification of malondialdehyde (MDA), glutathione (GSH), superoxide dismutase (SOD) and catalase (CAT) in midbrain tissue. Rotenone (ROT) exposure induced a significant increase (**p < 0.01) in the MDA level (a) and a decrease in the GSH level (b) in the midbrain of the ROTchallenged rats compared to the vehicle injected control (CONT) rats. NRD supplementation to the ROT-administered rats significantly decreased ( # p < 0.05) the level of MDA, and increased ( ## p < 0.01) the level of GSH. ROT injection also significantly decreased (**p < 0.01) the activity of SOD (c) and CAT (d) compared to the CONT rats. NRD supplementation significantly inhibited ( # p < 0.05) the ROT-induced decrease in SOD and CAT activity compared to the ROT-injected rats. Values are expressed as the mean ± SEM (n = 6-8) dopaminergic neurodegeneration, oxidative stress, inflammation, and α-synucleinopathy [6,22,24,25]. In the present study, the ROT-induced PD model was used to investigate the neuroprotective potential of the NRD phytochemical, which recently showed potential antioxidant activity [15]. ROT at a dose of 2.5 mg/kg of body weight once daily was used for 4 weeks to induce dopaminergic neurodegeneration. Our results show the novel potential benefits of NRD in attenuating dopaminergic neurodegeneration, improving antioxidant enzymes, and inhibiting inflammatory mediators and lipid peroxidation in the brain.
Tyrosine hydroxylase (TH) is a rate-limiting enzyme that catalyzes the synthesis of the DA neurotransmitter, which is synthesized in the dopaminergic neurons of the SNc area. DA is stored in synaptic vesicles and, in response to stimuli, released in the striatum to exert its physiological function [26]. In PD patients, more than 70 % of the dopaminergic neurons are dead resulting in the retraction of dopaminergic nerve terminals in the striatum and depletion of DA in the striatum [27]. Therefore, the death of dopaminergic neurons in the SNc is believed to decrease the level of striatal dopamine. The loss of dopaminergic neurons in the SNc and reduced density of striatal nerve terminals are considered to be the pathological hallmark of human PD [28].
In the present study, immunohistochemical examination of TH revealed that ROT causes a significant loss of TH+ dopaminergic neurons in the SNc area and the density of dopaminergic nerve terminal fibers in the striatum, which concurs with previous reports [22,25]. Interestingly, NRD supplementation to ROT-treated rats significantly protected against the ROT-induced loss of dopaminergic neurons in the SNc and reduced striatal nerve fiber density in the striatum. This result indicates that NRD has neuroprotective effects that prevent dopaminergic neuron loss  (green, a) and Iba-1-positive microglia (green, b) was found in the ROT-administered rats compared to the vehicle-injected CONT rats. In contrast, NRD supplementation to the ROT injected rats showed moderate staining of GFAP and Iba-1 compared to the ROT-injected rats (scale bar = 200 µm). Quantitative analysis of activated astrocytes (c) and microglia (d) revealed a significant increase (**p < 0.01) in the number of activated astrocytes and microglia in the ROT group rats compared to the CONT group. However, NRD supplementation significantly reduced ( # p < 0.05) the number of activated astrocytes and microglia in the ROT + NRD group compared to the ROT group. CONT rats and NRD-only-injected rats did not show any marked difference in the activation of astrocytes and microglia. Values are expressed as the percent mean ± SEM (n = 3)

Fig. 5
Proinflammatory cytokines in the midbrain tissue of the CONT, ROT, ROT + NRD, and NRD-only groups. The level of IL-1β (a), IL-6 (b), and TNF-α (c) was significantly increased (**p < 0.01) in the ROT-treated group when compared to the CONT group. However, NRD treatment significantly decreased ( # p < 0.05) the ROT-induced increase of these proinflammatory cytokines in the ROT + NRD group. No significant differences were observed in these cytokines between the CONT and NRD-only group. Values are expressed as the mean ± SEM (n = 6-7) and striatal nerve terminal retraction in ROT-induced PD in rats. A convincing number of studies have shown that dopaminergic neurons exist in a state of constant oxidative stress, as brain cells have a low antioxidant defense capacity and a tendency to generate reactive oxygen species due to the presence of highly oxidizable DA [29,30]. In our study, NRD administration protected against ROT-induced dopaminergic neurodegeneration possibly through its strong antioxidant action, as NRD shows efficient antioxidant activity in the hippocampus [15].
Defects in complex-I of the mitochondrial electron transport chain leads to a massive release of free radicals, and consequently, cellular death [31]. ROT, which is an inhibitor of complex-I of the mitochondrial electron transport chain, upon systemic administration to rats, causes the loss of ATP production and subsequently generates reactive oxygen species [30]. Furthermore, reactive oxygen species induce the oxidation of polyunsaturated fatty acid in a process called lipid peroxidation, which is characterized by the formation of MDA, a major product of lipid peroxidation that forms adducts with DNA bases and proteins, thereby causing cellular damage. The increased level of oxidative damage to DNA, proteins, and lipids, and decreased level of antioxidants have been reported in the brain of PD patients [4]. Significantly increased levels of MDA in the ROT-treated rats were observed compared to the vehicle-injected control rats. Interestingly, NRD supplementation to the ROT-injected rats significantly lowered the levels of MDA compared to the ROT-injected rats. Notably, NRD ameliorates against oxidative stress due to lipid peroxidation in the mouse hippocampus [15]. The inhibition of lipid peroxidation by NRD in the current study may have been mediated through the detoxification of peroxy radicals and reactive oxygen species, which further supports an antioxidative role for NRD.
GSH is an important antioxidant and has a crucial role in scavenging hydrogen peroxide. The decreased levels of GSH in the brain might indicate a state of oxidative stress. Decreased levels of GSH have been observed in surviving dopaminergic neurons of the SNc of PD patients compared to age-matched controls [32]. In the present study, the GSH level was significantly lower in ROT-treated rats than in control rats. Decreased levels of GSH lead to oxidative damage to DNA, protein, and lipids in PD [32][33][34]. Administration of NRD significantly increased GSH levels compared to the ROT injected and vehicle-treated rats (Group II). The effect of NRD on GSH may involve the direct antioxidant effects of NRD or the prevention of ROT-induced GSH oxidation.
SOD and CAT are the main antioxidant enzymes involved in detoxifying free radicals. The loss of SOD activity further contributes to an increase in oxidative stress in ROT-treated rats. Complex-I of the mitochondrial respiratory chain is the major source of superoxide radicals through inhibition of the electron transport chain [28]. The diminished activity of SOD would be detrimental in the scenario when superoxide radical production is increased. Our present data showed decreased activity of SOD in ROT-challenged rats compared to control rats. The administration of NRD was beneficial, because this treatment efficiently restored the SOD activity in ROTinjected rats. This result corroborated a previous report showing an increase in SOD activity after NRD administration [15]. CAT is an enzyme that is responsible for catalyzing the decomposition of H 2 O 2 . The maintenance of a normal metabolism of reactive oxygen species is important for proper cell function in different body parts [35,36]. In the present study, CAT activity was significantly diminished in ROT-injected rats, while the administration of NRD significantly reversed the decreased activity of CAT in NRD-supplemented rats. NRD supplementation reduced oxidative damage, as observed by the decrease in lipid peroxidation, restoration of the GSH level, and activities of the antioxidant enzymes (SOD and CAT) following ROT administration in rats.
Furthermore, ROT-treated animals showed an increase in neuroinflammation that is triggered and sustained through various mechanisms. Dysfunction of complex-I of the mitochondrial respiratory chain due to ROT treatment increased the release of reactive oxygen species, which can trigger the activation of glial cells [37]. A large number of studies report the involvement of neuroinflammation in PD, which is largely characterized by an accumulation of activated microglia [38]. The activation of astrocytes and microglia results in increased expression of GFAP and Iba-1, respectively. Microglial activation also causes NO· overproduction via iNOS induction [39]. Robust microgliosis and increased expression of iNOS have also been reported in postmortem brain samples of PD patients [40,41]. Proinflammatory cytokines such as IL-1β, IL-6, and TNF-α, and enzymes such as iNOS and COX-2 appear to be involved in DA toxicity, and have been observed in the cerebrospinal fluid and postmortem brain tissue samples of PD patients [40,[42][43][44][45].
In addition, inflammatory processes that are associated with the increased expression of COX-2 and iNOS are involved in the cascade of deleterious events that leads to neurodegeneration in PD [46]. In line with this report, here we showed an increase in activated astrocytes and microglia with a concomitantly enhanced level of the proinflammatory cytokines IL-1β, IL-6, and TNF-α in ROT-treated rats, which was significantly attenuated following NRD supplementation. Additionally, NRD supplementation normalized the expression of COX-2 and iNOS in the ROT-injected rats. Therefore, the present findings further suggest that NRD exerts its anti-inflammatory effects through inhibiting the increased level of COX-2 and iNOS, as well as proinflammatory cytokines. The NRD-mediated suppression of inflammation, as observed in the present study, could be facilitated by the strong anti-inflammatory activity of NRD that has been reported in other studies [47,48]. NRD has also been shown to inhibit the lipopolysaccharide-induced release of pro-inflammatory cytokines in bone marrow-derived dendritic cells, and diminish the induction of inflammatory mediators in an in vivo acute model of inflammation [49,50]. Oxidative stress and neuroinflammation play major roles in different neurodegenerative diseases such as stroke, Alzheimer's disease, and PD [51]. Drugs with anti-oxidative and anti-inflammatory characteristics are considered effective therapeutic agents for these neurodegenerative diseases [52]. Therefore, NRD, with its effective antioxidant and anti-inflammatory properties, may be responsible for the neuroprotection in the ROTinduced rat model of PD.

Conclusions
In conclusion, our study suggests that NRD supplementation ameliorates the ROT-induced dopaminergic neurodegeneration through its antioxidant and antiinflammatory activities. The present results highlight the therapeutic potential of NRD in neurodegenerative diseases including PD. However, further study is required to fully understand the molecular mechanisms of how NRD potentiates antioxidant and anti-inflammatory activities.


extfieldsofstudy: Biology Medicine

text:
Nuclear locus divergence at the early stages of speciation in the Orchard Oriole complex

Abstract As two lineages diverge from one another, mitochondrial DNA should evolve fixed differences more rapidly than nuclear DNA due to its smaller effective population size and faster mutation rate. As a consequence, molecular systematists have focused on the criteria of reciprocal monophyly in mitochondrial DNA for delimiting species boundaries. However, mitochondrial gene trees do not necessarily reflect the evolutionary history of the taxa in question, and even mitochondrial loci are not expected to be reciprocally monophyletic when the speciation event happened very recently. The goal of this study was to examine mitochondrial paraphyly within the Orchard Oriole complex, which is composed of Icterus spurius (Orchard Oriole) and Icterus fuertesi (Fuertes' Oriole). We increased the geographic sampling, added four nuclear loci, and used a range of population genetic and coalescent methods to examine the divergence between the taxa. With increased taxon sampling, we found evidence of clear structure between the taxa for mitochondrial DNA. However, nuclear loci showed little evidence of population structure, indicating a very recent divergence between I. spurius and I. fuertesi. Another goal was to examine the genetic variation within each taxon to look for evidence of a past founder event within the I. fuertesi lineage. Based on the high amounts of genetic variation for all nuclear loci, we found no evidence of such an event – thus, we found no support for the possible founding of I. fuertesi through a change in migratory behavior, followed by peripheral isolates speciation. Our results demonstrate that these two taxa are in the earliest stages of speciation, at a point when they have fixed differences in plumage color that are not reflected in monophyly of the mitochondrial or nuclear DNA markers in this study. This very recent divergence makes them ideal for continued studies of species boundaries and the earliest stages of speciation.

Introduction
For the last 25 years, neutral molecular markers, especially mitochondrial DNA (mtDNA), have been used to examine closely related specieselucidating their evolutionary histories and looking for evidence of gene flow. Studies of closely related taxa often address species delimitation by looking for reciprocal monophyly in mitochondrial genes (Moritz 1994;Zink and McKitrick 1995;Funk and Omland 2003;Hebert et al. 2004;Zink and Barrowclough 2008). Mitochondrial DNA has two key advantages that make it useful in this regard: (1) mtDNA has a rapid mutation rate, allowing for increased chances of accumulating variation; (2) mtDNA also has faster sorting (genetic drift) due to its maternal inheritance and generally lower effective population sizes (ESSs; Avise 1994;Zink and Barrowclough 2008).
However, there are many examples where species delimitation based solely on mtDNA would disagree with boundaries delimited using other characters (Baker et al. 2003;Funk and Omland 2003;Olsson et al. 2005;Omland et al. 2006;Knowles and Carstens 2007;Joseph and Omland 2009). There are several possible reasons for this disagreement. For instance, individual gene trees of closely related taxa can have unpredictable amounts of variation due to stochasticity in mutation, genetic drift, and sampling error (Rosenberg and Nordborg 2002;Knowles and Carstens 2007). A major factor that limits the utility of mtDNA is that because it is haploid and maternally inherited, different mitochondrial genes are inherited as a group, sharing a single evolutionary history that may or may not accurately represent the evolutionary history of the species of interest (Knowles and Carstens 2007). A second problem with the approach mentioned above is its focus on one criterionthe presence of reciprocal monophyly. Recent studies have shown that a substantial time lag exists between species divergence and the development of reciprocal monophyly in gene trees (Hey 2006;Omland et al. 2006;Knowles and Carstens 2007). Neutral genetic changes can diverge slower than traits under selective pressure, which may result in genetically indistinct taxa with fixed differences in phenotypes that play a role in fitness or sexual selection (e.g., Joseph and Omland 2009;Ross 2014). Species delimited by genetic markers thus can disagree with those delimited using other characteristics, such as morphology. Thus, recently diverged species will most likely not have any fixed differences in their mtDNA due to this time lag, as is evidenced by the multiple examples of known gene tree paraphyly in Animalia (see the following papers for a detailed explanation and for examples: Edwards and Beerli 2000;Funk and Omland 2003;Joseph and Omland 2009) and the many well-recognized species that lack fixed mtDNA differences (Ross 2014). Our study focuses on a pair of taxa that represent a well-documented case of mitochondrial paraphyly: the Orchard Oriole complex.
The Orchard Oriole complex is comprised of two taxa: Icterus spurius (Orchard Oriole), which breeds from southeastern Canada, across the eastern United States and into north-central Mexico, and Icterus fuertesi (Fuertes' Oriole), which breeds in a geographically restricted region along the coastal lowlands of southern Tamaulipas and Veracruz, Mexico (Howell and Webb 1995). These two taxa are considered separate species by the IOC World Bird List (Gill and Donsker 2015;also Navarro and Peterson 2004). Clements World Checklist (e.g., Clements 2007) formerly considered them separate species, but they were recently reclassified as subspecies within I. spurius, which is also followed by the American Ornithologist's Union (Chesser et al. 2015). In this article, we follow the taxonomy of IOC (Gill and Donsker 2015). The controversy regarding the separation of Orchard Oriole and Fuertes' Oriole into different species is based on a variety of studies carried out over the last 10 years. The two taxa differ in a number of ways. Adult males show a fixed difference in color (Fig. 1), with no overlap in color variation when their plumage is analyzed through spectrophotometric comparisons Kiere et al. 2007). Adult males of the two taxa are easily distinguished visually in the field (pers. obs.). They also differ in migratory behavior, as I. spurius migrates whereas I. fuertesi is a short-distance migrant or partial migrant (Jaramillo and Burke 1999;Tob on-Sampedro and Rojas-Soto 2015). These taxa breed at different latitudes, which has resulted in bioclimatic niche differentiation for their breeding distributions (Martin and Omland 2011). Lastly, studies of their vocalizations have found no evidence of song differentiation (Hagemeyer et al. 2012), yet a call that appears to play a role in territory defense differs significantly in a number of different measures of frequency, duration, and amplitude (Sturge et al. 2016).
In addition, Baker et al. (2003) sampled these two taxa across their ranges and sequenced two regions of mtDNA (cytochrome b and control region). They found that the two species were paraphyletic in their mtDNA with intermixed haplotypes (haplotypes that are shared and/or show paraphyletic relationships with another species; Joseph and Omland 2009). Icterus spurius and I. fuertesi show both these patterns of intermixing in their mtDNA. However, the frequencies of different mitochondrial haplotypes differ significantly, indicating little or no mitochondrial gene flow between the two taxa (F ST = 0.608, P < 0.00001, Baker et al. 2003). Baker et al. (2003) suggested that retained ancestral polymorphisms could be sufficient to explain this mitochondrial paraphyly between the two taxa, but they did not test it with rigorous statistical analyses (i.e., coalescent methods using multiple loci).
The mitochondrial paraphyly in the Orchard Oriole complex has two possible explanations, as with other cases. First, the split between the two taxa could have happened so recently that the haplotype variation within each of these taxa has not yet coalesced; thus, they could share haplotypes due to retained ancestral polymorphisms. Second, it is possible that gene flow has been occurring since the two taxa diverged, which could explain the observed patterns of haplotype and lineage sharing. In avian species, hybridization and gene introgression are common even among more distantly related species Price 2008;Rheindt and Edwards 2011), and species divergence can occur even in the face of ongoing gene flow (Hey 2006).
Additionally, the results from Baker et al. (2003) suggested that the divergence that led to the formation of these two taxa could have been the result of a founder event within the lineage leading to I. fuertesi. This was based on two mitochondrial genes (cytochrome b and control region), for which I. fuertesi showed lower haplotype diversity than I. spuriuswhich is the predicted pattern if the lineage leading to I. fuertesi was founded from a few individuals whose genetic diversity only represented a portion of the available diversity within the ancestral population (Harrison 1991;Baker et al. 2003;Funk and Omland 2003). The two taxa within the Orchard Oriole complex exhibit different migratory behavior, as mentioned abovewith I. spurius migrating long distances while I. fuertesi is a short-distance migrant that overwinters in the southern portion of its breeding range, where it likely overlaps with overwintering I. spurius (Howell and Webb 1995;Jaramillo and Burke 1999;Kondo and Omland 2007; Tob on-Sampedro and Rojas-Soto 2015). As migratory behavior is highly labile in orioles (Winker 2000;Kondo and Omland 2007), a loss or decrease in migratory behavior within a subset of the ancestral population could have resulted in the founding of the current breeding grounds of I. fuertesi, followed by peripheral isolates speciation in allopatry (Mayr 1942(Mayr , 1963West-Eberhard 2005). Examining Figure 2, I. fuertesi looks like it could be a classic example of peripheral isolates speciation based on its much smaller breeding range, which is geographically isolated from the closest I. spurius breeding populations by the Sierra Madre Oriental Mountains. However, the lack of genetic variation within I. fuertesi reported by Baker et al. (2003) could also have resulted from differences in sampling effortthere were roughly twice as many I. spurius individuals as I. fuertesi included in their study.
We addressed these questions regarding the mitochondrial paraphyly between these taxa, and the potential founder event involving I. fuertesi by sequencing multiple loci (four nuclear introns) and through an increase in sample size, adding nine I. fuertesi samples. We sequenced multiple nuclear loci to obtain a broader picture of the genomic divergence between these two taxa. Using multiple loci, we attempted to avoided either stochastic or deterministic processes that might have only affected a part of the genome or one type of marker (e.g., mtDNA; Edwards and Beerli 2000;Peters et al. 2008). We also used population coalescent methods, including IM (isolation with migration, Hey and Nielsen 2004;Hey 2006) to examine polyphyly and shared haplotypes to see whether it was possible to differentiate between retained ancestral polymorphisms and gene flow as possible explanations for these patterns. Our goals were to address the following questions: (1) Do nuclear loci also support very recent divergence, or does mtDNA reflect recent introgression between two well-diverged lineages? (2) What is the nature of the nuclear intron variation between these two taxa? Do they support two distinct groups? (3) Was divergence the result of a founder event and peripheral isolates speciation? (4) Is it possible to determine whether these two taxa are currently exchanging genes?

DNA sequencing
In this study, we had 32 samples already available to us and we collected nine additional I. fuertesi samples from multiple localities in Veracruz, Mexico, to increase our sample size. Not all of the old samples would amplify; as a result, we included samples from 25 I. spurius and 14 I. fuertesi from across their breeding ranges (Fig. 2, Table S1). For each individual, we sequenced a total of five loci: four nuclear introns and one mitochondrial gene. The four nuclear introns were located on four different chromosomes of the zebra finch genome; three were autosomal loci (GADPH11 397 b.p., TGFB2 574 b.p. and RDP2 298 b.p.) and one was on the Z sex chromosome (SLC9 429 b.p.). The fifth locus was the mitochondrial control region (CRdomain I; 344 b.p.). Accession numbers on GenBank for the genes we sequenced are as follows: For additional information about the five loci, see Table 1.
We extracted DNA from muscle tissue of all samples using a DNeasy Blood and Tissue Kit (QIAGEN, Valencia, CA). PCRs were then performed on each locus using a Gene Amp PCR System 9700 (Applied Biosystems, Foster City, CA), and we verified our product using 1% agarose gels with ethidium bromide. We cleaned PCR product using QIAquick PCR Purification Kits (QIAGEN). The amplified loci were sequenced using ABI's BigDye v. System 2400 (Applied Biosystems). We used EDTA-ethanol precipitation to remove excess dye terminators. Chromatograms for each sequence were produced using the University of Maryland, Baltimore County's ABI3100 DNA Sequencer (Applied Biosystems). We edited and aligned the chromatograms using Sequencher 4.1 (Gene Codes Corporation, Ann Arbor, MI).
Individuals within our dataset that contained more than one polymorphic site within a locus were phased using PHASE 2.1.1 (Stephens et al. 2001;Stephens and Scheet 2005). We set the burn-in rate to 1000 and the thinning interval to 1. Each run included 10,000 iterations, and we repeated each analysis ten times with different random numbers as the starting seed. We selected the haplotype states suggested by the runs with the highest log-likelihood values and included all haplotypes that resolved at a probability of at least 0.5 (for justification of this approach and cutoff value, see Jacobsen and Omland 2012). All of the SNP phasing in this study met this cutoff, so no further action was necessary. TOPALi 2.5 (Milne et al. 2004) was used to look for evidence of intralocus recombinationnone of our loci had a significant DSS peak and so all were included in our dataset. We constructed and edited haplotype networks within NETWORK v. 4.6.1.1 (Bandelt et al. 1999) using the median joining algorithm.

Population statistics and structure
The number of haplotypes (h), the haplotype diversity (Hd), the nucleotide diversity (p), the average number of nucleotide differences per site (p), and the number of All the dots for I. fuertesi represent localities where we collected multiple samples. This map also shows the breeding ranges for both species. The map was built using ArcView 3.2 (Environmental Systems Research Institute, Redlands, CA) and the distributional ranges from Nature Serve (www.natureserve.org) and by adding each locality point where the samples were collected. Chromosome locations were determined using BLAST searches of a reference genome (Zebra Finch). The number of alleles for each taxon are given as N.
nucleotide polymorphisms (h) were all calculated using DnaSP (Rozas et al. 2003; for results, see Table 2). We also used DnaSP to calculate Tajima's D (Tajima 1989), Fu and Li's F* (Fu and Li 1993), Ramos-Onsins & Roza's neutrality state (R 2 , Ramos-Onsins and Rozas 2002), and mismatch distributions (with associated Tau and Theta values). To look for significant population structure within the nuclear loci, we used PGDSpider to produce a dataset that could be analyzed within Structure by assigning each haplotype a unique number (Lischer and Excoffier 2012). One of the assumptions of Structure is that the different loci included are unlinked. As mentioned previously, our loci were located on different chromosomes of the Zebra Finch genome, so our approach should not violate this assumption. Note that this method does not consider the similarity between different haplotypes when it assigns numerical designations. We analyzed our data using Structure v. 2.3.4 (Pritchard et al. 2000;Hubisz et al. 2009). We used the admixture model with correlated allele frequencies and the LOCPRIOR flag activated, and set both the Dirichlet parameter for degree of admixture (a) and the allelic frequency parameter (k) to be inferred from the dataset, rather than being fixed. All runs had a burn-in of 50,000, followed by 250,000 iterations to collect data to estimate posterior probabilities. We tested K values from 1 to 10, and replicated each run 20 times for all K values. We used Arelquin v3 3.5.1.2 (Excoffier and Lischer 2010) to run AMOVAs on each locus to test for significant population structure within and between the two taxa (Schneider et al. 2000;Excoffier and Lischer 2010). We calculated Φ ST values in Arelquin as well, which are analogous to Wright's F ST values, but take into account the number of mutations between each haplotype within the sample (Excoffier and Lischer 2010).

Coalescent analyses
To determine whether or not we could estimate how much of the shared DNA variation between these two taxa is due to retained ancestral polymorphisms versus gene flow (either historic or ongoing), we analyzed mtDNA, nuclear DNA, and the combined mtDNA and nuclear DNA datasets using both the programs IM and IMa (isolation with migration, Hey 2006Hey , 2007. IM estimates seven demographic parameters: ESSs of the ancestral (h A ) and daughter populations (h 1 and h 2 ), migration rates between the two daughter populations (m 1 and m 2 ), and the time the two daughter populations diverged from each other (t), and the contribution of the ancestral population to each daughter population (splitting parameter, s) (Hey and Nielsen 2004). We used the infinite sites model and set the inheritance scalars based on the modes of inheritance for the loci (autosomal were set to 1, Z linked to 0.75 and mitochondrial to 0.25). In these analyses, we included multiple heated chains (30 chains) and examined the ESS estimates to insure the Markov chains were mixing appropriately (Hey and Nielsen 2004). We also included 1,000,000 generations and sampled every 100 generations so we could monitor output. We completed the run multiple times, using different starting seeds for each run. We also varied the following parameters in an effort to reach convergence: the maximum migration rates for both populations (m 1 and m 2 ), the maximum time of population splitting (t), the duration of the burn (b) and the duration of the run (l), the heating parameters (g 1 and g 2 ), the number of chains (n) and the number of chain swap attempts per step (k), and the estimates of maximum ancestral and daughter population sizes (q A , q 1 , and q 2 ).

Results
Both taxa were polymorphic with respect to all loci included in this study, as evidenced by the allele networks (Fig. 3). None of the nuclear introns show species monophyly, and there are many shared alleles between taxa. In contrast, the mitochondrial CR is much more segregated between the taxa, with only one shared haplotype and very little intermixinghowever, two I. spurius haplotypes are only one base pair different than the closest I. fuertesi haplotype, while they are either two or three base pairs different than the closest I. spurius haplotypes.
Both taxa show high levels of haplotype diversity (Hd) for all of the nuclear loci, ranging from 0.724 (SLC9 for I. spurius) to 0.889 (TGFB2 for I. fuertesi) ( Table 2). However, both show lower levels of diversity for the mitochondrial gene CR, with 0.597 for I. spurius and 0.242 for I. fuertesi. Note that I. spurius has more than two times the haplotype diversity of I. fuertesi for CR, yet the taxa show similar amounts of diversity for all four nuclear loci (Table 2). In terms of nucleotide diversity (p), three of the five loci (CR, TGFB2, and RDP2) have higher nucleotide diversity for I. spurius and two (GADPH11 and SLC9) have higher diversity for I. fuertesi. The majority of the loci for both taxa show negative values for Tajima's D and for Fu and Li's F*, with a few exceptions (Table 3). However, none of these values differed significantly from zero once Bonferroni corrections were applied, indicating that there is no evidence of selection within the loci included in this study. The sizes of the Ramos-Onsins & Rozas R 2 values indicate that both taxa may have undergone recent population expansions (Ramos-Onsins and Rozas 2002).
The results from STRUCTURE analysis of the four nuclear loci provide no evidence of population structure. The highest log-likelihood probability for Ln P(D) was detected for K = 1 (Fig. 4), and structure was unable to group individuals into clusters that matched taxon boundaries (Fig. 5). The AMOVA results for these loci generally agree with this finding, as they show that the majority of the haplotype variation can be explained by within population divergence (see Table 4). Thus, there is generally a lack of nuclear intron structure between the two taxa, with a percent of variance explained by species boundaries that ranges from À6.3% (SLC9) to 5.6% (RDP2), see Table 4 for additional values. These findings are further supported by three of the four nuclear loci, which had nonsignificant Φ ST values, with the vast majority of the variation being explained by within population variation. The remaining autosomal locus, TGFB2, however, has an important difference; whereas the majority of variation within this locus is again explained by within population variation, TGFB2's Φ ST value was significant, with a value of 0.056 (P < 0.00001).
For the IM analyses of our nuclear loci, the IM runs consistently failed to converge even after more than twenty runs with the autosomal locithe last five of which included the Z-linked intron SLC9and using multiple random number seeds. (The failure to converge could be show extensive haplotype sharing between Icterus spurius and I. fuertesi. In contrast, the haplotype network for the mitochondrial control region shows very little intermixing and only one shared haplotype between the two taxa. Circle sizes correspond to sample sizes of each allele. caused by recent divergence; other recent publications from our laboratory used similar methods and resulted in convergence; see Discussion). As a result, we will not present the results from our IM analyses in this article. The results of the AMOVA for mtDNA CR are a clear contrast to the nuclear DNA (Table 4). These results show statistically significant structure between the two taxa, with 71.3% of genetic variation being explained by species boundaries. Additionally, the Φ ST value for this locus is 0.717 (P < 0.00001), indicating there is evidence of a great deal of structure between the two taxa for this locus. With the addition of more I. fuertesi, the mtDNA haplotype frequencies are clearly statistically significantly different.

Discussion
Do nuclear DNA support very recent divergence?
The majority of the nuclear loci included in this study show very little evidence of structure between the two taxa. One nuclear locus, TGFB2, is the exception as it has a small, but significant Φ ST value (Φ ST = 0.056, P < 0.00001). The amount of variation explained at the among-species level in the AMOVA for TGFB2 is 5.1%, indicating that a small portion of the genetic variation can be explained by species boundaries. If the mtDNA paraphyly first described by Baker et al. (2003) was the result of a recent mitochondrial introgression between two well-diverged lineages, we would expect to find a great deal of divergence between the nuDNA loci that would not be reflected in the mtDNA (Jacobsen and Omland 2012). Instead, these two lineages appear to be very early in the process of divergingas evidenced by the lack of structure for the majority of nuclear loci, as well as by the lack of reciprocal monophyly in mtDNA.
With increased sampling, we found significant differences in mitochondrial haplotype frequencies between the two populations that confirms gene flow between the two taxa may be restricted (Baker et al. 2003). For the mitochondrial gene CR, there is clear evidence of significant population structure between the two taxa, with a Φ ST value of 0.717 (P < 0.00001). Even with the inclusion of eight more I. fuertesi samples, only one haplotype is shared between the two taxa for this locus, and only two haplotypes are intermixed and showing mtDNA paraphyly. Two haplotypes that we found only in I. spurius are only one base pair different from I. fuertesi haplotypes (Fig. 3, CR), yet are at least two base pairs different than the closest I. spurius haplotypes. Thus, with increased sampling effort, CR seems to be nearing reciprocal monophyly between groups. These lineages show a classic pattern of intermediate divergence in which the only haplotype shared between the taxa is a central, likely ancestral haplotype. Omland et al. (2006) discussed this as an early step on the road to monophyly ("neotypy") indicating that these two taxa are clearly showing evidence of divergence within the faster sorting mtDNA.
Does nuclear DNA support two distinct groups?
In comparison, the nuclear DNA included in this study show very little evidence of structure between the two taxa. The IM analyses failed to converge even though we tried a range of strategies to produce convergence. We   have published a number of recent studies from our laboratory on other species of recently diverged Icterus using similar methods, which resulted in convergence in those cases (Jacobsen and Omland 2012;Cort es-Rodr ıguez et al. 2013;Cort es-Rodr ıguez and Omland 2016). Furthermore, in the Orchard species complex, all four nuclear loci show a great deal of intermixing of alleles (Fig. 3). Additionally, the results of the AMOVAs indicates that almost all of the variation is found within the populations of each taxon (Table 4), further supporting a lack of divergence within the slower sorting nuclear loci included in the study. This is what we would expect to find if the divergence between these two taxa was at the very earliest stages of speciation.
Although the majority of variation within the nuclear loci exists at the level of variation within populations, one autosomal locus showed a small amount of variation that could be explained by species boundaries. TGFB2 had roughly 5% of its variation being explained by variation within each taxon (P < 0.00001). Thus, one of the nuclear loci shows some evidence of segregation, indicating that the nuclear genomes of the two taxa are beginning to diverge.
Our inability to find much structure for the nuclear loci does not necessarily mean the taxa are not evolutionarily discrete units. There are many examples of known mitochondrial paraphyly within ornithological literature, with a large portion of it being attributed to retained ancestral polymorphisms (Funk and Omland 2003;Mckay and Zink 2010). Neutral nuclear markers, with their larger ESSs and slower sorting rates, are predicted to lag even further behind (Price 2008;Zink and Barrowclough 2008;Joseph and Omland 2009). Taxa that fall into this category may still be discrete in other ways that indicate they are evolutionarily distinct units, for example, for traits under sexual selection such as song or plumage coloration (Edwards et al. 2005;Marthinsen et al. 2008;Price 2008;Joseph and Omland 2009;Mckay and Zink 2010). A study examining a recent divergence in two species of crows, Corvus (corone) corone and C. (corone) cornix, showed significant divergence in the genes that controlled for plumage coloration of the two taxa, yet a lack of divergence in neutral genetic markers due to hybridization and widespread introgression, showing that assortative mating and sexual selection can lead to genetic differentiation in regions under selection (Polestra et al. 2014). Similar to the above study, research in our laboratory has previously shown that I. spurius and I. fuertesi have fixed differences in adult male coloration Kiere et al. 2007). In a parallel study to this one, we found fixed differences between the taxa in a call that seems to play a role in territory defense, thus might be under the influence of sexual selection (Sturge et al. 2016). Any genes that could be contributing to either plumage coloration or this call thus have the potential to have diverged between these two taxa.
Did I. fuertesi result from a founder event?
If the divergence event that separated I. fuertesi and I. spurius resulted from a few individuals of the migratory common ancestor founding a new, nonmigratory or short-distance migratory population, followed by peripheral isolates speciation (Mayr 1942(Mayr , 1963West-Eberhard 2005), we would predict that I. fuertesi would have significantly less genetic diversity than I. spurius. Looking at both the haplotype networks (Fig. 3) and the population statistics for each locus (Table 2), it is clear that both I. spurius and I. fuertesi have similar levels of genetic diversity for all nuclear loci included in this study (although the mitochondrial locus still indicates less genetic diversity within I. fuertesi). Based on these results, we found no evidence of a founder event in I. fuertesi's recent past. This result is surprising because the distribution of I. fuertesi suggested that it might have evolved through this mechanism. Instead, both taxa show a great deal of allelic diversity, indicating relatively large historic population sizes in spite of the very restricted range of I. fuertesi. Nevertheless, it is conceivable that I. fuertesi was established through a large number of founders, and/ or that subsequent gene flow from I. spurius has added a lot of genetic diversity through introgression. As the IM analyses did not converge, we are unable to completely eliminate this possibility without including many additional loci. However, a classic splitting of a single more widespread population clearly fits our data. This taxon pair would be ideally suited to next generation sequencing approaches to study this and other aspects of their divergence (e.g., testing for genes affecting plumage coloration and ecological tolerances). Whole genome comparisons would most likely provide more insight into both the divergence between these two taxa, and any potential gene flow that has occurred since this divergence. The statistical results of the Ramos-Onsins and Rozas R 2 for both taxa (Table 3) does seem to support recent and rapid population expansion for both taxa, which may have contributed to problems with IM (Hey and Nielsen 2007). Nevertheless, both taxa also have a great deal of both nucleotide and haplotype diversity, as well as unique haplotypes, supporting the existence of large historic ESSs.
Is it possible to determine whether the two taxa are currently exchanging genes?
Based on the lack of convergence in IM, due largely to the lack of structure within the nuclear loci, we are unable to determine how much of the shared variation is the result of retained ancestral polymorphisms or is the result of gene flow. The control region shows far more evidence of structure, with only one shared haplotype and two intermixed haplotypes within all individuals included in this study. There are two possible explanations: (1) significantly different haplotype frequencies in mtDNA between the taxa indicates gene flow is likely greatly reduced or nonexistent and is not being captured at the nuDNA level due to retained ancestral polymorphisms, or (2) haplotype frequency differences for only the mtDNA could also be explained by Haldane's rule in the face of ongoing gene flow (Mckay and Zink 2010;Peters et al. 2012). Haldane's rule predicts that hybrids of the heterogametic sex will have reduced fitness or be sterileso that mtDNA in birds would be less likely to introgress from one species to another (Coyne 1985;Peters et al. 2012). As a consequence of this, as well as the failure of the IM analyses to converge, it is difficult if not impossible to determine whether gene flow is currently occurring between these two taxa, or has historically occurred since their divergence. Although their populations breed in allopatry, gene flow could occur as migratory I. spurius pass through I. fuertesi's breeding range in early spring. However, the significant population structure, not only within the mtDNA (which is less likely to introgress), but also within the nuclear locus TGFB2, indicates that I. spurius and I. fuertesi are in the process of diverging genetically from one another. Based on this and other evidence such as coloration and migratory patterns, it is likely that these two taxa have become evolutionarily distinct groups, regardless of any potential gene flow.

Conclusions -Model Example of Early Divergence
Icterus spurius and I. fuertesi are two taxa that are in the very earliest stages of speciation. Their mtDNA show only a few shared and intermixed haplotypes, having almost achieved reciprocal monophyly for this locus. The four nuclear loci included in this study show very little population structure between the two taxayet one locus (TGFB2) shows a small, but significant Φ ST divergence between the two taxa, confirming the formation of two distinct genome pools. Based on these molecular findings, as well as the fixed differences in plumage, migratory behavior, and bioclimatic niches, there is strong evidence that these two taxa are just beginning to diverge from one another, making them ideal study organisms for the early stages of speciation and the formation of species boundaries Kiere et al. 2007;Kondo and Omland 2007;Martin and Omland 2011).


s2fieldsofstudy:
Biology
Engineering
Materials Science

text:
Patterning Biomaterials for the Spatiotemporal Delivery of Bioactive Molecules

The aim of tissue engineering is to promote the repair of functional tissues. For decades, the combined use of biomaterials, growth factors (GFs), and stem cells has been the base of several regeneration strategies. Among these, biomimicry emerged as a robust strategy to efficiently address this clinical challenge. Biomimetic materials, able to recapitulate the composition and architecture of the extracellular matrix, are the materials of choice, for their biocompatibility and higher rate of efficacy. In addition, it has become increasingly clear that restoring the complex biochemical environment of the target tissue is crucial for its regeneration. Toward this aim, the combination of scaffolds and GFs is required. The advent of nanotechnology significantly impacted the field of tissue engineering by providing new ways to reproduce the complex spatial and temporal biochemical patterns of tissues. This review will present the most recent approaches to finely control the spatiotemporal release of bioactive molecules for various tissue engineering applications.

FiGURe 1 | Schematic showing a spatially patterned fibrous material functionalized with different sets of delivery systems (A), in separate compartments (B) (Minardi et al., 2014). Schematic showing the temporal patterning of a material with two sets of delivery systems (C), for the staged release of bioactive molecules (D). engineering applications (Balint et al., 2014), including bone and cartilage regeneration (Tampieri et al., 2008;Henkel et al., 2013;Minardi et al., 2015). The biomaterials market for implantable devices is estimated to be worth $33 billion by the end of 2019 (Transparency Market Research, 2012). Tissue healing in musculoskeletal surgery has benefited from the surgical implantation of bio-conductive scaffolds (Hsu et al., 2012). However, the use of GFs and chemo/cytokines in clinical practice yielded controversial results (Garner et al., 2011;Anitua et al., 2012). One famous example involved the use of recombinant human bone morphogenetic protein-2 (rhBMP-2) in spinal fusion (Medtronic, INFUSE ® : rhBMP-2-infused collagen scaffold). Although approved by the US Food and Drug Administration (FDA) in 2002, the high doses of rhBMP-2 used by this device resulted in major adverse side effects, including cancer, spinal cord compression from soft tissue swelling, spinal cord impingement from ectopic bone formation, elevated bone resorption from osteoclast activation, and preferential induction of adipogenesis over osteogenesis (Kaneko et al., 2000;Smucker et al., 2006;Wong et al., 2008;Robin et al., 2010;Zara et al., 2011;Epstein, 2013).
These side effects were caused by the massive release of rhBMP-2 from the collagen sponge after implantation (5-10 mg of BMP-2/implant) (McKay et al., 2007). This dose is over 5000 times higher than the amount required for bone formation in vitro (Wang et al., 1988;Place et al., 2009) and even higher than physiologic doses found in the osteogenic niche where other osteogenic stimuli (chemical, physical, and structural) additionally contribute to induce bone formation (Place et al., 2009;Tampieri et al., 2011). When the uncontrolled pharmacokinetics showed adverse side effects, the FDA revoked the device's approval (McKie et al., 2014).
After this, GF dose-response effects became all the more crucial (Shields et al., 2006), and it became evident that adverse outcomes could have been avoided with the controlled, localized release of rhBMP-2. It was determined that burst release and widespread tissue exposure to GFs, together with an early dissipation from the implantation site (e.g., with irrigation, bleeding, and edema) were not ideal pharmacokinetics for tissue regeneration (Lai et al., 2013). Furthermore, it has become clear that the first-order release of a single bioactive molecule is not sufficient to mimic the complex biochemical gradients present at a specific stage of tissue regeneration (Wang et al., 2013;Minardi et al., 2014).
Extraordinary progress has been made toward the design of biomaterials with suitable multiscale hierarchical structures, facilitating the staged release of a combination of bioactive molecules, according to any complex delivery pattern (Biondi et al., 2008;Guldberg, 2009;Chen et al., 2010).
Bioactive factors can be incorporated within materials using layer deposition or integrated into their fibrous mesh via electrospinning or self-assembly techniques (Sun et al., 2003;Hosseinkhani et al., 2006;Minardi et al., 2014). Herein, we review the most recent strategies for biomaterial fabrication, featuring spatiotemporal patterns of bioactive molecules.

enGineeRinG 3D BiOPATTeRneD MATeRiALS
The three-dimensional (3D) interactions between cells and the extracellular matrix (ECM) have been proven to be crucial to orchestrate tissue formation and regeneration in response to injury. Materials such as hydrogels and scaffolds engineered to emulate the ECM can support tissue healing (Hoffman, 2012;Loh and Choong, 2013).
Three-dimensional patterning has been defined as the entrapment of biochemical (Lee et al., 2015a) or structural heterogeneity (Zorlutuna et al., 2011) within the structure of a material and successfully employed in the design of biomimetic materials for tissue engineering. To emulate tissues' biochemical gradients, the patterning of bioactive molecules within biomaterials should be controlled over time and space (Figure 1). This has been accomplished using multiple strategies that facilitate controlled spatial and temporal release kinetics.
Hydrogels capable of promoting cell viability and interactions with elements of the ECM have been advantageous for many regenerative medicine applications (Nguyen and West, 2002). Hydrogels are highly hygroscopic polymers that can be fully engineered to mimic specific chemical and physical properties of tissues and are therefore often used as scaffolding materials for tissue engineering applications (Tibbitt and Anseth, 2009). Additionally, due to their biocompatibility and ability to control the release rates of bioactive molecules, hydrogels have been successfully used as drug reservoirs in tissue engineering applications (Peppas et al., 2006). The 3D patterning of hydrogels is crucial to enhance their spatial heterogeneity and improve their features for cell seeding. 3D patterning can fall into two categories: those that integrate patterning during fabrication (i.e., stereolithography) or those that involve post-processing of a uniform hydrogel (Ahmed, 2015). Horn et al. (2007) created hydrogels suitable for spinal cord repair by combining thiol-modified hydroxyl-apatite with acrylate-functionalized poly(ethylene glycol) (PEG). Tsang et al. (2007) developed a multilayer PEG hexagonal hydrogel functionalized with hepatocyte cells and demonstrated better cell viability in the center of the construct compared to uniform hydrogel disks. Lim and Sun (1980) used alginate-based polymers cross-linked with calcium ions for the treatment of diabetic animals. De Souza et al. (2009) demonstrated the in vivo biocompatibility of chitosan, phospholipids, and lauric aldehyde or lauric chloride hydrogel blends. Hydrogels can also be patterned by indirect binding. In an elegant recent study, Stupp and coworkers described a successful example of this strategy: they proposed peptide amphiphile nanofibers gels with binding affinity for BMP-2 (Rajangam et al., 2008;Lee et al., 2015b), to solve the health concerns associated with the use of this factor in patients. Others have successfully functionalized the surface of gels and scaffolds with heparin or others glycosaminoglycans to dock bioactive molecules in situ (Liang and Kiick, 2014;Corradetti et al., 2016).
One of the major concerns in the engineering of materials functionalized with bioactive molecules for regenerative medicine is the preservation of the payload . Due to their structure, 3D scaffolds proved successful in minimizing bioactive factors' exposure to harsh conditions in vivo, preserving their bioactivity (Yuan and Liu, 2012).
Scaffolds have been biochemically patterned through many different approaches. An established strategy consists of the integration of microparticles into the biomaterial structure (Chen et al., 2010). Integrating delivery systems within the scaffold matrix enables the design of scaffolds with patterns of various geometries and purpose, while controlling their release (López-Noriega et al., 2015).
Similarly, our group proposed a chitosan-gelatin scaffold functionalized with composite microspheres consisting of mesoporous silicon microparticles and poly(lactic-co-glycolic acid) (PLGA) for the controlled release of small bioactive molecules (Pandolfi et al., 2016). Wei et al. (2007) showed that incorporating PLGA nanospheres into nanofibrous poly(l-lactic acid) (PLLA) scaffolds loaded with rhBMP-7 enhanced osteogenesis. Kim et al. (2004) accomplished the incorporation of a combination of drugs into a PLGA-based fibrous mat, fabricated through electrospinning.
Interface tissue engineering focuses on the development of tissue grafts capable of replacing defective interfaces, such as ligament-to-bone, tendon-to-bone, and cartilage-to-bone (Almodóvar et al., 2014). These interfaces exhibit anisotropic structural properties, which gradually vary from one tissue to another. Using homogeneous biomaterials ultimately leads to graft failures (Seidi et al., 2011). Singh et al. (2008) introduced an interesting microparticle-based scaffold fabrication technique as a method to create 3D scaffolds with spatial control over multiple bioactive molecules using uniform PLGA microspheres. They demonstrated that embedding the PLGA microparticles into their scaffolds led to more sustained payload release. Recently, we proposed a multiscale approach to selectively integrate different types of nanostructured composite microspheres in a multicompartmented collagen scaffold (Minardi et al., 2014). By fully embedding the microspheres in the type I collagen matrix of the scaffold, the authors were able to spatially pattern the microspheres into the different compartments, and the collagen coating on the microspheres allowed for the zero-order release of the payload for almost 2 months. This method of scaffold functionalization proved to preserve the bioactivity of cytokines for several weeks. The controlled release of the payloads over long periods of time favors on-scaffold cell recruitment while avoiding adverse effects due to the bolus administration of therapeutic molecules in the tissues surrounding the implant . To achieve the combined release of multiple molecules from different compartments, Song et al. (2012) incorporated silica nanoparticles into electrospun fibers. Detamore's group attempted osteochondral regeneration through the creation of continuous gradients of two bioactive factors within a 3D scaffold (Dormer et al., 2010). The authors developed a polymer microsphere-based scaffold, which could be fabricated using microspheres loaded with various GFs to create continuous and opposite gradients of BMP-2 and transforming growth factorbeta (TGF-β) (Mohan et al., 2011).

CeLL ReSPOnSe TO 3D BiOCHeMiCALLY PATTeRneD MATeRiALS
Most of the current understanding of the regenerative process is based on the use of simplified in vitro 2D cell cultures, which fail to reproduce tissue complexity. Recently, 3D cultures on biomimetic materials mimicking the ECM microenvironment are becoming the in vitro models of choice to study the regenerative process of specific tissues. 3D biomimetic cultures are more physiologically relevant than those in 2D and simpler than in vivo models (Hoffmann and West, 2010). As reviewed in the previous paragraph, over the past decade various technologies have been developed to create spatiotemporal gradients for complex biomaterials (Santos et al., 2012). Graded materials can elicit various cell behaviors, such as adhesion, orientation, motility, surface antigen display, cytoskeletal condensation, activation of tyrosine kinases, and modulation of intracellular signaling pathways that regulate transcriptional activity and gene expression (Figure 2). Cells are able to decode information provided by the ECM and respond to specific stimuli, such as topography (Ishii et al., 2005;Raghunathan et al., 2013), mechanical properties (Ulrich et al., 2009;Humphrey et al., 2014), bioactive signals (Li and Folch, 2005;Yañez-Soto et al., 2013;Malik et al., 2015), and concentration gradients of both soluble and tethered GFs (Gattazzo et al., 2014). An established strategy to elicit such biological functions consists of the functionalization of hydrogels with biochemical cues, as described above (Tibbitt and Anseth, 2009;Geckil et al., 2010). Synthetic biomaterials, such as hydrogels or bulk polymeric scaffolds, demonstrated low cellular adhesion, necessitating functionalization with bioactive molecules (Tian et al., 2012). Since synthetic hydrogels lack biochemical cues, different strategies for their 3D biochemical patterning have been developed to increase their level of biomimicry (Luo and Shoichet, 2004) and enhance their interaction with cells (Yan et al., 2011). To favor cell adhesion, migration, and differentiation within hydrogels, the main strategy has been to immobilize adhesive peptides (e.g., RGD) (Bellis, 2011). The use of short peptides to tune cellular response to hydrogels showed success due to the stability of the chemical conjugation between the peptide of interest and the material, without altering its conformation.
Recently, in a very elegant study, the Anseth research group proposed an RGD patterned hydrogel, synthesized through a versatile process of sequential bio-orthogonal click-chemistry reactions (DeForest et al., 2009). Their approach was further developed by García and coworkers, who proposed light to trigger adhesive peptides on the surface of their hydrogel in vivo (Lee et al., 2015c). Adhesion and proliferation are not the only two important features to achieve tissue homeostasis, but also cell migration is fundamental for all morphogenetic processes in tissue regeneration (Sternlicht and Werb, 2001). Cell migration occurs in response to gradients of soluble or insoluble signals. In tridimensional constructs, migration is more difficult to achieve due to the mechanical resistance of the surrounding ECM. Synthetic materials must contain cell-adhesive ligands for traction, GFs as signals of migration, as well as space to allow cell movement. Taraballi and colleagues synthesized a plethora of different hydrogels functionalized with bioactive peptides (BMPH1 andBMPH2) that allowed for the adhesion, proliferation, and differentiation of neural stem cells into neurons both in vitro and in vivo (Gelain et al., 2010(Gelain et al., , 2011Taraballi et al., 2010). They demonstrated that neural cells were able to migrate inside the gel but only extended neurites when both factors were presented at the same time. Their findings suggested that several signaling molecules could work best in combination, benefiting from their synergistic effect.
However, in the natural ECM, cells not only respond to signals presented on the surface of the ECM but also to soluble stimuli, especially in the regeneration process. Usually, these factors diffuse through the ECM and bind to their specific receptors on cells' surface, activating specific transduction cascades. In vivo, the transient nature of GF signaling is combined with the slow, sustained signals received from the ECM (Lund et al., 2009). For example, epidermal growth factor (EGF), FGFs, TGF-β, and platelet-derived growth factor (PDGF) (Farokhi et al., 2013) have been used to accelerate wound healing by inducing both epithelial cell and fibroblast proliferation, as well as de novo matrix deposition (Kim et al., 2012). Similarly, the release of insulin-like growth factor 1 (IGF-1) and TGF-β1 from polymer scaffolds functionalized with delivery systems showed to successfully induce chondrogenic differentiation (Ertan et al., 2013).

COnCLUSiOn
In the past decade, significant advances have been accomplished in the design of biopatterned materials able to accomplish temporally and/or spatially controlled release of bioactive molecules. Despite these advances, methods should be further developed to prepare patterns and gradients with controlled shape and kinetics, in order to tune the desired cell mechanisms in vivo. Altogether, the studies herein reviewed show the potential of 3D biomaterials spatiotemporally patterned with bioactive molecules to recapitulate the complex biochemical milieu of target tissues.

AUTHOR COnTRiBUTiOnS
SM and FT conceived the idea of this review. SM wrote the Abstract, Introduction, and Conclusion sections. LP drafted the first paragraph, FT drafted the second paragraph, and SM and FT revised and finalized the mini-review with the help of ET.

ACKnOwLeDGMenTS
Authors thank Ms. Megan Livingston for helping with editing the text.

FUnDinG
This work was supported by the Brown Foundation (Project ID, 18130011), the Hearst Foundation (Project ID, 18130017), by the Cullen Trust for Health Care Foundation (Project ID, 18130014) and by DoD USAMRMC (Project ID, W81XWH-15-1-0718).