s2fieldsofstudy: Biology Chemistry

text:
Enhanced autophagy contributes to excitotoxic lesions in a rat model of preterm brain injury

Cystic periventricular leukomalacia is commonly diagnosed in premature infants, resulting from severe hypoxic-ischemic white matter injury, and also involving some grey matter damage. Very few is known concerning the cell death pathways involved in these types of premature cerebral lesions. Excitotoxicity is a predominant mechanism of hypoxic-ischemic injury in the developing brain. Concomitantly, it has been recently shown that autophagy could be enhanced in excitotoxic conditions switching this physiological intracellular degradation system to a deleterious process. We here investigated the role of autophagy in a validated rodent model of preterm excitotoxic brain damage mimicking in some aspects cystic periventricular leukomalacia. An excitotoxic lesion affecting periventricular white and grey matter was induced by injecting ibotenate, a glutamate analogue, in the subcortical white matter (subcingulum area) of five-day old rat pups. Ibotenate enhanced autophagy in rat brain dying neurons at 24 h as shown by increased presence of autophagosomes (increased LC3-II and LC3-positive dots) and enhanced autophagic degradation (SQSTM1 reduction and increased number and size of lysosomes (LAMP1- and CATHEPSIN B-positive vesicles)). Co-injection of the pharmacological autophagy inhibitor 3-methyladenine prevented not only autophagy induction but also CASPASE-3 activation and calpain-dependent cleavage of SPECTRIN 24 h after the insult, thus providing a strong reduction of the long term brain injury (16 days after ibotenate injection) including lateral ventricle dilatation, decreases in cerebral tissue volume and in subcortical white matter thickness. The autophagy-dependent neuroprotective effect of 3-methyladenine was confirmed in primary cortical neuronal cultures using not only pharmacological but also genetic autophagy inhibition of the ibotenate-induced autophagy. Strategies inhibiting autophagy could then represent a promising neuroprotective approach in the context of severe preterm brain injuries.

Introduction
The important progress done in neonatal care constantly increases the survival rates of premature infants. Conversely, the proportion of neurological disabilities developed by survivors is hardly reduced especially for those with severe impairment. One of them is the diplegic cerebral palsy, called also spastic diplegia of Little 1 , affecting still between 3-7% of very low birth weight (VLBW) infants 2,3 . The strongest predictor of this form of cerebral palsy in VLBW infants is cystic periventricular leukomalacia (cPVL) 4 , a form of preterm white matter (WM) injury adjacent to the lateral ventricles, which occurs either from a hypoxic-ischemic (HI) event around birth or after infectious events such as septic shock, necrotizing enterocolitis or even reported after viral infections 5,6 . Improving the outcomes for these severely affected babies remains a challenging health issue.
Beside inflammation and reactive oxygen species formation, excitotoxicity seems to be crucial in the pathophysiology of many preterm brain injuries such as PVL 3,7,8 . Excitotoxicity consists in an excessive or prolonged activation of excitatory amino acid receptors (especially those of glutamate) due to a failure of sufficient reuptake and/or excessive release at the synaptic level. It induces a massive increase in intracellular calcium concentration and thus activates numerous intracellular cascades potentially leading to neuronal cell death 9,10 . Glutamate homeostasis is highly important for human brain development (proliferation, migration, differentiation, survival processes and synapses refinement) 11 . However it also confers to immature brain a vulnerability to excitotoxic injuries since a higher level of ionotropic glutamate receptors are expressed in developing brain compared to that of adult 8,[12][13][14][15] . These receptors are in addition more readily activated. Excitotoxic lesions can occur following a panel of deleterious events (that can be related) such as infection/inflammation, hypoxia and/or ischemia. Excitotoxicity is then a common pathological mechanism of various perinatal brain injuries. In neurons the peak of expression of NMDA receptors appears to occur at term in which grey matter (GM) damage is predominant than in preterm 13 . In human WM, this peak occurs in preterm brain glial cells, especially in pre-oligodendrocytes O4 + , and contributes to the high sensitivity of preterm WM. PVL was mainly thought to be associated to WM injury but it is clearly shown now that GM damage is also often involved in a "neuronalaxonal disease" 16,17 .
Experimental research has revealed the complexity of the pathophysiology of excitotoxic death showing multiple and interrelating cell death mechanisms reflected by mixed features of neuronal death including not only the well-known "apoptotic-necrotic continuum" with features of apoptosis and necrosis in the same dying neurons 18 but also characteristics of enhanced macroautophagy 10,[19][20][21] . Autophagy is a physiological cellular mechanism of degradation and recycling of dysfunctioning long lived proteins and organelles 22 . Its main form (macroautophagy, hereafter called autophagy), consists in the formation of a multimembrane intermediate compartment, named autophagosome, that engulfs part of the cytosol containing proteins and organelles to be degraded. The autophagosome then fuses with a lysosome, forming an autolysosome, to degrade its content through lysosomal hydrolases activity 22 . Autophagy is thus essential for cellular homeostasis and could be used as a survival response to different stresses such as nutrients deprivation, accumulation of toxic proteins or pathogen invasion 23 . However, dysregulated increase in autophagic process has been also implicated in cell death as an independent mechanism (termed "autophagic cell death") or more frequently as a mediator of other types of cell death, mainly apoptosis, and then designed as "autophagymediated cell death" 10,[24][25][26][27][28] . Abnormal high level of autophagosomes and autolysosomes with increased lysosomal enzyme activity were often observed in dying neurons in models of excitotoxicity including perinatal cerebral HI 19,20,[29][30][31] . Interestingly, we also recently demonstrated excessive autophagic features in postmortem brains of human term newborns presenting severe hypoxic-ischemic encephalopathy (HIE) 29,31 . Although controversies remain concerning the role of autophagic activation 24,[32][33][34][35] , most of the studies using autophagy inhibition, either through pharmacological inhibitors such as 3-methyladenine (3-MA) 30,[36][37][38][39][40] or through specific and genetic inhibition of autophagyrelated genes (atg) 20,29,31,41 , have revealed a pro-death role of autophagy in perinatal and adult cerebral HI models.
The present study aims to determine the role of autophagy in excitotoxic lesions of the premature brain using a widely recognized rodent model that mimics some features of cPVL 42 . Autophagy flux and the neuroprotective effect of autophagy inhibition either pharmacologically with 3-MA or genetically in neuronal cultures was investigated in the context of an excitotoxic insult induced by an injection of the glutamate analogue ibotenate. An involvement of autophagy in excitotoxic preterm brain damage would reveal a new cellular death pathway and open the way for new neuroprotective strategies in severe preterm brain injuries.

Results
Ibotenate injection in subcortical WM of rat pups induces brain injury As a model mimicking some aspects of cPVL, we selected a widely recognized model of preterm excitotoxic brain injury consisting in applying an intracerebral injection of ibotenate in rat pups 42 . The injection of ibotenate (10 µg) in the subcortical WM at the level of the right cingulum of 5-days-old rat pups causes a severe brain damage of the WM and GM resulting in ventricular enlargement as illustrated by cresyl-violet-stained coronal sections in Fig. 1a 24 h and 16 days after the insult. Immunoblot analysis of two death markers, the calpaindependent 150-145 kDa spectrin/fodrin (SPTAN) fragment and cleaved CASPASE-3 (CASP3), indicated an activation of both calcium-dependent necrotic cell death and caspase-dependent apoptosis (Fig. 1b) at 6 and 24 h after the insult (~5 and 14 fold increases respectively at 24 h). Although the profile of CASP3 activation was similar between females and males, we observed 24 h after the insult a greater activation (~2 times more) in females. Double immunolabeling showed that CASP3 was activated in neurons (RBFOX3-positive cells) 24 h after ibotenate injection (Fig. 1c).
To evaluate the effect of ibotenate on WM, the thickness of the subcortical WM was measured at three different levels around the ibotenate injection site. Reduced thickness could reflect an indirect effect of the GM lesion with secondary axonal degeneration and/or a direct excitotoxic effect on immature oligodendrocytes. A significant decrease in the thickness of subcortical WM of 13% for corpus callosum, 40% for the WM located at the level of cingulum and 20% for external capsule was detected in ibotenate-treated brains (Fig. 5d). These WM alterations were not observed when 3-MA was co-injected with ibotenate (Fig. 5d).
Taken together, these data showed a long term beneficial effect of pharmacological autophagy inhibition against ibotenate-induced brain damage that was moreover observed both in males and females.

An excitotoxic dose of ibotenate enhances autophagic flux in primary cortical neuronal cultures
In order to confirm that ibotenate could enhance autophagy and that specific autophagy inhibition could be neuroprotective against ibotenate-mediated excitotoxicity, we then decided to evaluate the effect of an excitotoxic dose of ibotenate on autophagy in primary cortical neuronal cultures. In general, 50 µM of ibotenate promoted neuronal death as suggested by a 3 and 14 fold increase in lactate deshydrogenase (LDH) release in the culture medium 3 and 6 h, respectively after ibotenate treatment (Fig. 6a). This LDH release was abolished by MK801 or EGTA pretreatment (data not shown) confirming excitotoxic mechanisms (NMDA receptors overactivation and calcium overload).
As shown in Fig. 6b, activation of autophagy occurred along with neuronal death since both LC3-II expression and SQSTM1 degradation were increased at 3 and 6 h. In order to confirm an enhanced autophagic flux, ibotenate was first applied in the presence of lysosomal enzymes inhibitors (Fig. 6c). A combination of E64 and pepstatin A1 (PepA) (10 μg/ml) induced an accumulation of both LC3-II (of 123%) and SQSTM1 (of 13%) reflecting the failure in the autophagy degradation step. When ibotenate was applied 4 h after E64/PepA, LC3-II and SQSTM1 accumulations were even greater (174 and 21% respectively), demonstrating that ibotenate treatment triggered the new formation of autophagosomes and thus confirming that ibotenate treatment increased the autophagosome biogenesis. Second, autophagic flux was monitored using the tandem mRFP-GFP-LC3 plasmid that allows to discriminate between LC3 expressed in neutral compartments (GFP + RFP+; early autophagosomes: i.e. autophagosomes before fusion with lysosomes) and in acidic vesicles (GFP-RFP+; i.e., late autophagosomes: autophagosomes after fusion with lysosomes (autolysosomes)) thanks to the pH sensitivity differences exhibited by the two fluorescent proteins (Fig. 6e). A quantification of the different LC3-positive dots per neuron per µm 2 clearly demonstrated that both autophagosomes formation (~5 fold increase of GFP + RFP + dots) and their fusion with the lysosomes (~13 fold increase of GFP−RFP + dots) were enhanced by ibotenate treatment. Taken together, these results on cortical neuronal cultures allow to conclude that an excitotoxic dose of ibotenate was efficient to induce a boost of neuronal autophagic flux.

Pharmacological and genetic inhibition of autophagy is protective against ibotenate-induced excitotoxicity in primary cortical neuronal cultures
We then assessed the functional role of the ibotenateenhanced autophagy in primary cortical neuronal cultures. We first used 3-MA. Pre-treatment with 3-MA (10 mM) prevented the increases in both LC3-II expression and SQSTM1 degradation at 6 h after ibotenate treatment (Fig. 7a). 3-MA displayed a significant neuroprotective effect as shown by a decrease of 26% in LDH release (Fig. 7b). Interestingly, blocking lysosomal degradation with E64/PepA also reduced neuronal death (Fig. 7c) suggesting a role of autophagy degradation step in ibotenate-induced neurotoxicity.
Then, in order to inhibit more specifically autophagy, downregulation of the expression of two important autophagy proteins (ATG7 and BECLIN1 (BECN1)) was performed using lentiviral vectors transducing short hairpin RNAs (shRNAs) (Fig. 8a). Transduction of cultured primary cortical neurons with Atg7 and Becn1 shRNAs resulted in an efficient inhibition of autophagy as demonstrated by a decrease in both LC3-II expression and SQSTM1 degradation (Figs. 8b, c). The death-promoting role of enhanced autophagy in ibotenate-induced neuronal death was confirmed by a reduction of~30% of LDH release when ATG7 and BECN1 were downregulated (Fig. 8d).
In conclusion, these in vitro data confirmed that an excitotoxic dose of ibotenate could induce autophagymediated neuronal death and that targeting autophagy inhibition could lead to neuroprotection.

Discussion
In the present study, we investigated for the first time the role of autophagy in a preclinical model of premature brain damage related to cPVL 42 . The effect of ibotenate injection on rodent brain development is highly   Values are mean ± SEM, *p < 0.05; **p < 0.01; ***p < 0.001 dependent of the age. If injected around P5, when neuronal migration is completed, the model mimics some of the preterm brain injury features 42 . Human preterm WM injury is often observed with GM developmental alteration and/or damage (neuronal-axonal disease) 2,3,16,17,44 . Persistent cerebral volume reduction and ventriculomegaly are also observed in premature infants compared to full-term 45,46 . The injection of ibotenate in the subcortical WM at the level of the right cingulum caused a severe lesion in rat pups and led to reduced subcortical WM thickness, significant loss of tissue and important lateral ventricle enlargement 16 days after injury. It has been shown that excitotoxicity-induced inflammation through microglial release of cytokines and free radicals played a central deleterious role 47,48 and that apoptotic pathways were implicated in similar models of preterm brain injury 49,50 . However, autophagy-mediated neuronal death has never been investigated in this context.
Macroautophagy is an important physiological mechanism of degradation present at a basal level complementary to the proteasome system. A "controlled" upregulation of autophagy has been considered for a long time as a survival response, for instance, acting as an alternative source of energy during starvation or as a quality control step eliminating toxic metabolites, defective organelles or intracellular pathogens 23,51,52 . However, this dogma has been challenged about 30 years ago by the description of dying cells without typical morphological hallmarks of apoptosis or necrosis and containing numerous autolysosomes 27 . This new morphological type of cell death was coined "autophagic cell death" and the missuse of the term led to a strong debate on a possible pro-death role of autophagy in some stress 8.6 ± 1.4%; Ibo + 3-MA: 74.7 ± 7.5%). c E64/PepA pretreatment reduces neuronal death as indicated by a decrease in LDH release measured 6 h after Ibo addition (veh DMSO: 16.8 ± 2.9; Ibo DMSO: 100 ± 5%; veh E64/PepA: 21.3 ± 5.8%; Ibo E64/PepA: 69.6 ± 9.7%). Values are mean ± SEM, *p < 0.05; **p < 0.01; ***p < 0.001 conditions 10,34,35 . In fact, "autophagic cell death" is a distinct mechanism of cell death independent of apoptotic and necrotic machinery that is observed in some specific circumstances 27,38,53 . However, there is now evidences and it is well accepted, that autophagy can be more frequently involved in cell death as a trigger leading to   20,21,29,34,38,[54][55][56][57] . Excitotoxicity and cerebral HI are some of the conditions where autophagy is enhanced 58 . Although some controversies exist 32,33,59 , most of the studies using pharmacological inhibitors such as 3-MA 30,37 and especially those using specific and genetic inhibition of atg 20,29,31 revealed a deleterious role of autophagy. They also supported the concept of a strong interconnection between autophagic and apoptotic mechanisms in perinatal cerebral HI 19,20,32,60 . The present study is the first to demonstrate a death-promoting role of autophagy in a preterm model of excitotoxic brain lesion.
Our results strongly suggest that autophagy is enhanced in neurons after ibotenate injection in the brain of rat pups, as also demonstrated in primary cortical neurons cultures. Autophagosome formation is increased from 6 h as shown by a higher level of LC3-II and more LC3positive dots. This increase was not due to impaired autophagosome degradation that would have occurred if lysosomal function was defective. In fact, p62/SQSTM1 is also degraded and autolysosomes are increased in neurons as indicated by larger and more numerous LAMP1-and CTSB-positive dots. This result is consistent with EM observations showing the presence of both autophagosome-like multimembrane vesicles and autolysosome-like dense structures in dying neurons of ibotenate-injected rat pups. We provide here, especially in primary neuronal cultures, compiling evidences demonstrating that the autophagy flux is increased by ibotenate treatment. Both Western blots against LC3 and p62/ SQSTM1 with or without E64/pepstatinA co-treatments and the use of the GFP-RFP-LC3 construct lead to the same conclusion, ibotenate treatment is increasing autophagosome formation and degradation. Moreover, the GFP-RFP-LC3 construct we used is less sensitive to acidification comparing to other constructs, such as the mKate2-pHluorin-hLC3 61 , suggesting that the number of autolysosomes could be theoretically underestimated. Our in vitro data on primary cortical neuronal cultures also strongly argued for the induction of an autophagy-mediated neuronal death by ibotenate. Inhibition of autophagy, not only pharmacologically with 3-MA or E64/PepA, but also genetically by downregulating two important ATG proteins, ATG7 or BECN1, was neuroprotective as previously shown when excitotoxicity was induced in hypoxic conditions 29 .
Even if caution is necessary concerning the limited specificity of 3-MA as an autophagy inhibitor 22 , we here clearly showed that the dose used in vivo efficiently prevented autophagy (decrease in LC3-II level and SQSTM1/p62 degradation). Furthermore, since it is known that permanent impaired autophagy leads to neurodegeneration (as it would be the case with long term genetic inhibition of autophagy) 43 , the use of 3-MA was appropriate to study long term effect on brain lesion in vivo (16 days after the insult). When ibotenateenhanced autophagy was prevented by 3-MA treatment, both CASP3 activation and calpain-dependent cleavage of SPTAN (used as an indicator of necrotic cell death characterized by a calcium increase) were significantly decreased. Moreover, ultrastructural observations showed condensed chromatin in nuclei of dying and highly autophagic neurons. Mixed features of apoptosis and autophagy were also observed in rodent models of perinatal cerebral HI 19,20,30,31,60 and we previously demonstrated that autophagy could contribute to apoptosis using widely recognized apoptotic stimuli in primary cortical neuron cultures 62 . The beneficial effect of 3-MA on neuronal cell death, including apoptosis, resulted here in a strong significant neuroprotective effect at long term on both GM and WM, suggesting a crucial role of autophagy in mediating (apoptotic and necrotic) cell death. Studies on perinatal cerebral HI models have also shown that autophagy inhibition could be neuroprotective and reduce apoptosis 20,30,31 suggesting that the transient inhibition of autophagy could be a promising strategy to protect the immature brain against excitotoxic insults.
It has been shown that females are more resistant to perinatal cerebral injury than males, especially in the context of cerebral HI and after moderate lesions, in rodents as well as in humans [63][64][65][66] . The reason of this gender difference is still unclear, but sex-dependent cell death pathways have been recognized after perinatal cerebral HI, especially more active caspase-dependent pathways in females [66][67][68] . In the present study, we also found a stronger CASP3 activation and a more important variability in females in almost all the different parameters investigated. However, mean ibotenate-induced brain lesion volume was similar in both genders. Moreover, the protective effect of 3-MA against both ibotenate-induced autophagy and cerebral lesions was as efficient in males as in females.
In conclusion, we showed for the first time that enhanced autophagy could mediate cell death in a premature model of excitotoxic brain damage. Autophagy inhibition in this severe model is very promising since the protective effect obtained is similar or even better to other previously described neuroprotectants such as caspases inhibitors 69 , erythropoietin 70,71 , BDNF 72,73 , melatonin 74 or magnesium sulfate 75,76 . Interestingly, we recently demonstrated that autophagy is enhanced in dying neurons in the ventrolateral nucleus of the thalamus and the lentiform nucleus of term newborns with severe HIE 31,60 . Apoptotic markers were also expressed in dying highly autophagic neurons, arguing for a possible association between autophagy and apoptosis also in humans. Moreover an increased presence of autophagosomes (LC3-positive dots) was recently demonstrated in WM injury of extremely preterm infants 77 . These 3 different studies using human newborn brain sections 31,60,77 and the important neuroprotection obtained with 3-MA in the present study in a severe model of cPVL suggest that enhanced neuronal autophagy could be a promising target. The development of strategies transiently inhibiting autophagy could pave the way for new therapies against neonatal severe excitotoxic injuries such as HIE and cPVL.

Primary cortical neuronal cultures
Primary neuronal cultures were prepared according to the Swiss laws for the protection of animals from pieces of cortices of 2-day-old Sprague-Dawley rat pups (Janvier Labs, Mayenne, France). The procedures were approved by the Vaud Cantonal Veterinary Office. After dissection, dissociation and trituration, neurons were plated in neurobasal medium (Gibco, NY, USA; 21103-049) supplemented with 2% B27 (Gibco; 17504044), 0.5 mM L-glutamine (Sigma, MO, USA; 49419) and 100 μg/ml penicillin-streptomycin (Gibco; 15140122) and maintained at 37°C with a 5% CO 2 -containing atmosphere as described previously 78 . Western blot analyses were done on neurons plated at a density of~7 × 10 5 cells/dish (35-mm poly-D-lysine pre-coated dishes (BD Biosciences, NJ, USA; 356467) and at a density of~3 × 10 5 cells on 12-mm glass coverslips coated with 0.01% poly-L-lysine (Sigma; P4832) for immunocytochemistry and imaging. For all the Western blots or imaging results, at least three independent experiments, each involving two or three culture dishes or coverslips, were performed.

Quantification of cell death with lactate dehydrogenase release
Cell death was assayed by measurement of LDH released in the medium using the Cytotox 96 nonradioactive cytotoxicity assay kit (Promega, WI, USA; G1780) as previously described 78 . LDH measurements were normalized with respect to the values of ibotenatetreated neurons 6 h after the ibotenate addition.

Knockdown of ATG using lentiviral vectors
Downregulation of Atg genes were performed with pLKO lentiviral vectors (Open Biosystems/Dharmacon, CO, USA) expressing rat-specific shRNA sequences from TRC (the RNAi consortium) library as described previously 62 . A combination of TRCN0000092163 and TRCN0000092166 for Atg7 (GenBankTM NM_001012097), TRCN0000033552 for Becn1 (Gen-BankTM NM_053739.2) and a pLKO vector containing scrambled shRNA (Open Biosystems/Dharmacon) as control vector were used. Primary cortical neuron cultures were infected at DIV7 with 50 ng of the viral capsid protein p24/ml culture medium for each vector.

Rat model of preterm excitotoxic brain injury
All experiments were performed in accordance with the Swiss laws for the protection of animals and were approved by the Vaud Cantonal Veterinary Office. Ten µg of ibotenate (diluted 5 µg/µl in acetic acid 0.02%) were stereotaxically injected under isoflurane anesthesia (2.5%) in the subcortical WM at the level of the right cingulum (1 mm posterior and 1 mm right from Bregma and 1.5 mm depth from the skull surface) of 5-day-old male and female Sprague Dawley rats (Janvier Labs) (model adapted from Marret and colleagues (1995) 42 ). The control animals received an injection of the same volume of vehicle (acetic acid 0.02%). The pharmacological autophagy inhibitor 3-MA (2 µl of 30 mg/ml in saline) was stereotaxically injected in the right lateral ventricle (0.5 mm posterior and 1 mm right from Bregma and 2.5 mm depth from the skull surface) just before ibotenate injection. Control animals received an injection of the same volume of saline as vehicle. After recovering from anesthesia, rat pups returned to the dam until sacrifice.
For immunofluorescence labeling, Alexa Fluor donkeyanti-rabbit or mouse secondary antibodies (Invitrogen; A21202, A21203, A21206, A21207) were used. A LSM 710 Meta confocal microscope (Carl Zeiss) were used for confocal laser microscopy. Images were processed with LSM 510 software and mounted using Adobe Photoshop.

Quantification of autophagic and lysosomal labeling
Confocal images of fluorescent immunostaining against LC3, CTSB and LAMP1 were acquired using the LSM 710 Meta confocal laser scanning microscope (Carl Zeiss) and images were then processed with Adobe Photoshop CC 2015. Positive dots were quantified using ImageJ software and expressed as a number of positive dots per neuron per μm 2 and, for the lysosomal markers CTSB and LAMP1, as a mean dot area per neuron per μm 2 .

Cerebral regions volume and WM thickness measurements
Sixteen days after the injury (at P21) brains were perfused, frozen and entirely cut into series of 20 μm coronal sections spaced at 500μm disposed in series. On a cresyl violet-stained series, the ipsilateral volume of the total surviving tissue, the cortex and the lateral ventricle were measured using the Zen Blue software (Zeiss). The volumes were then expressed as a percentage of the total brain volume. WM thickness were measured on 3 consecutive cresyl violet-stained sections starting from the first on which the genu of the corpus callosum appeared (approx. 0.6-0.8 mm anterior to the Bregma according to the "atlas of the rat brain in stereotaxic coordinates at P21" of Khazipov et al., http://www.ialdevelopmentalneurobiology. com/images/atlases/Atlas-p21.pdf 81 ). The thickness of the ipsilateral corpus callosum (on the midline) and the subcortical WM at the level of the cingulum (1.4 mm apart from the midline) and at the beginning of the external capsule (2 mm apart from the midline) were measured parallel to the midline with the Zen Blue software (Zeiss). Values are expressed as a mean of the three measures.

Statistics
Values were expressed as mean values ± standard error of the mean (SEM). Data were analyzed statistically using GraphPad PRISM (version 7.03) software. The normality of the distribution was first tested using Shapiro-Wilk tests. Parametric data were analyzed using a Welch's ANOVA test (one-way ANOVA with unequal variances) followed by a post-hoc Tukey-Kramer test. Nonparametric data were analyzed using a Kruskal-Wallis test (non-parametric analog of the one-way ANOVA) followed by a post-hoc Steel-Dwass. P < 0.05 was chosen as threshold for statistical significance.


extfieldsofstudy: Medicine Biology

s2fieldsofstudy:
Biology
Environmental Science

text:
Genetic variability and ontogeny predict microbiome structure in a disease-challenged montane amphibian

Amphibian populations worldwide are at risk of extinction from infectious diseases, including chytridiomycosis caused by the fungal pathogen Batrachochytrium dendrobatidis (Bd). Amphibian cutaneous microbiomes interact with Bd and can confer protective benefits to the host. The composition of the microbiome itself is influenced by many environment- and host-related factors. However, little is known about the interacting effects of host population structure, genetic variation and developmental stage on microbiome composition and Bd prevalence across multiple sites. Here we explore these questions in Amietia hymenopus, a disease-affected frog in southern Africa. We use microsatellite genotyping and 16S amplicon sequencing to show that the microbiome associated with tadpole mouthparts is structured spatially, and is influenced by host genotype and developmental stage. We observed strong genetic structure in host populations based on rivers and geographic distances, but this did not correspond to spatial patterns in microbiome composition. These results indicate that demographic and host genetic factors affect microbiome composition within sites, but different factors are responsible for host population structure and microbiome structure at the between-site level. Our results help to elucidate complex within- and among- population drivers of microbiome structure in amphibian populations. That there is a genetic basis to microbiome composition in amphibians could help to inform amphibian conservation efforts against infectious diseases.

Introduction
Multicellular organisms act as hosts to a diverse suite of bacterial communities, collectively referred to as the microbiome. Recent research across multiple taxa has highlighted that these bacterial communities perform beneficial functions for the host, including protection from infectious pathogens [1][2][3]. Despite the importance of symbiotic relationships between microbes and their hosts, a comprehensive understanding of factors that influence the diversity and composition of the microbiome, particularly for non-human animals, is lacking. Though previous studies have revealed landscape-scale variation among populations in microbiome structure [4][5][6], and fine-scale variation among individuals within populations due to life stage [6], diet [7], gender [7] and genotype [7,8], few have investigated both within and among population drivers of microbiome structure in a single framework [6,8]. Studying factors that determine variation in the distribution of bacterial symbionts at both the individual and landscape scale is crucial for understanding how individual susceptibility to pathogens varies within and among populations.
The amphibian skin microbiome represents a model system for understanding the ecological drivers of host microbiome structure and the role of symbiotic bacteria in protecting the host from pathogens. The lethal and globally distributed amphibian chytrid fungus (Batrachochytrium dendrobatidis; Bd) is causing mass mortalities and rapid population declines of amphibians around the world and is a major driver of the current amphibian extinction crisis [9]. Certain host-associated bacteria have been shown to influence susceptibility to Bd, with particular bacterial genera conferring increased or decreased protection against Bd in experimental and field studies [2,[10][11][12]. Host genotype is one factor that may underpin variation in the presence and persistence of these protective microbes, both within and among host populations [7,8,[13][14][15][16]. For example, polymorphism of the major histocompatibility complex (MHC) Class IIb gene has been shown to co-vary with gut microbiome structure in vertebrates [7]. The expression of antimicrobial peptides (AMPs) on frog skin is governed by innate immunity [17], and so variation in immune genes may drive structural changes in the microbiome through differential AMP expression. The population genetic structure could therefore yield differences in microbiome at the landscape level if it is representative of the non-random distribution of functional genetic variation at immune loci. Furthermore, shifts in the environmental reservoir of bacteria across large geographic scales may also cause differences in the relative abundances of bacteria able to colonise the host [4]. Within populations, fine-scale variation in life-history traits, such as age, may further modify the ability of bacteria to coexist on the skin, or favour differential selection by the host from the environment. Though microbiome differences have been shown among larvae, juveniles and adults in some amphibian species [18,19], changes in microbiome structure throughout larval ontogeny have not been investigated. This may be especially pertinent, given that fine-scale variation in Bd infection loads has been found across developmental stages of tadpoles [20,21]. Research into the links between host genetics, developmental stage, microbiome composition and pathogen susceptibility may enable a better understanding of the factors governing disease in vulnerable populations.
Here we use a disease-challenged high-altitude amphibian species in the Drakensberg Mountains of southern Africa to investigate within and among population predictors of microbiome structure that may have implications for disease susceptibility. Using microsatellite genotyping and 16S amplicon sequencing, we examine how population genetic structure, individual genotype, and developmental stage determine differences in microbiome structure among individuals. Specifically, we aim to: (i) quantify Bd infection across Amietia hymenopus frog populations; (ii) characterise the microbiome of A. hymenopus; (iii) examine the genetic population structure of A. hymenopus on the Drakensberg plateau; (iv) investigate predictors of microbiome structure among sites (specifically river, geographic distance and genetic differentiation among host populations); (v) investigate within-site predictors of microbiome structure (specifically host genotype and ontogenetic stage); and (vi) examine the roles of host microbiome and host genetics in influencing Bd infection loads.

Study site
The Drakensberg Mountains of southern Africa are a vast mountain range reaching up to 3482 m above sea level (asl). Here, the Mont-aux-Sources plateau crosses north-east Lesotho and South Africa. On the South African side, four streams (Vemvane, Tugela, Bilanjil and Ribbon Falls) flow 1-3 km across the plateau from multiple origins in the high peaks on the South Africa-Lesotho border, before dropping off a~1000 m precipice into South Africa (Fig. 1). Two other streams on the plateau flow into Lesotho; Khubela and Kubedu. The Phofung river frog (Amietia hymenopus) is the only amphibian found in most of the Mont-aux-Sources stream system, with the exception of the Kubedu river, where A. hymenopus is absent and A. umbraculata is present [22]. Long-term (10 year) sampling in this region indicates A. hymenopus has a consistent history of Bd infection with associated mass mortalities [23].

Sample collection
We collected tadpoles from Mont-aux-Sources by permission of Ezemvelo KwaZulu Natal Wildlife and the Kingdom of Lesotho Ministry of Tourism, Environment and Culture. We sampled eight sites across five rivers over the plateau; 'Nampolice' and 'Vemvane' from the Vemvane River (South Africa); 'Tugela', 'Tugela 2' and 'Tukelahed' from the Tugela River (South Africa); 'Khubnam' from the Khubela River (Lesotho); and 'Bilanjil' and 'Ribbon Falls' from their respective rivers (South Africa; Fig. 1). Highresolution mapping was unavailable for the study area, so we used QGIS (version 2.18) to combine Open Source data sets and then digitise the water courses within the study area (Supplementary Information).
We collected 20-28 tadpoles per site (nine at Khubnam) along a~150 m transect within an altitudinal range of 3010-3090 m asl to minimize variation in Bd infection. We conducted two sampling trips 4-weeks apart; in March 2015, we collected tadpoles from six sites (Nampolice, Tugela, Tugela 2, Tukelahed, Bilajil and Ribbon Falls; Fig. 1) and analysed these for Bd infection, host genotype and microbiome composition. In April 2015, we collected tadpoles from Vemvane and Khubnam to enhance the population genetics analyses, but did not carry out microbiome analyses on these samples in order to avoid spurious variation caused by temporal changes (Table 1). We stored tadpoles in 95% ethanol prior to Gosner staging [24] and exported carcasses for molecular analyses in the UK by permission of the Department of Agriculture, Forestry and Fisheries, South Africa, and the Department for Environment, Food and Rural Affairs, UK. We excised full mouthparts for DNA extraction as due to their keratinised nature, Bd primarily affects this area in tadpoles. DNA was extracted using the DNeasy® Blood and Tissue Kit (Qiagen) and was used for all subsequent analyses.

Bd infection
We constructed serial dilutions of DNA and conducted quantitative PCRs on a Qiagen Rotor Gene according to Boyle et al. [25]. We quantified Bd infection intensity in triplicate using a 1:10 dilution of DNA with Bd standards ranging from 0.025 to 250 zoospores/µl. We calculated infection prevalence and intensity for each site.

16S amplicon sequencing of bacterial communities
We sequenced the microbial communities present on tadpole mouthparts using a modified version of the protocol in Kozich et al. and Bates et al. [26,27]. Transient bacteria most likely rinsed off during the ethanol preservation and so microbial analyses reflect the resident microbiome in the mouthparts of tadpoles. Briefly, we amplified a V4 region of the 16S rRNA gene in triplicate for 139 samples, as well as two negative controls (one per PCR plate). Triplicates were pooled and normalized using a titration run on an Illumina MiSeq using v2 nano chemistry, followed by a full run using 250 bp paired-end v2 chemistry. We also sequenced a mock community containing 20 bacteria isolated from amphibian skin (genera Acinetobacter, Citrobacter, Enterobacter, Flavobacterium, Plantibacter, Pseudomonas and Serratia).

Analyses
We performed sequence processing in DADA2 v1. 4 [31]. Following Longo and Zamudio [32], we excluded SVs that contained fewer than 100 reads across the entire dataset (n = 19,395), yielding 2384 SVs. Results using all SVs are similar and not shown here. After filtering, mean library size across individuals was 20,160 reads (range 13,156-83,063). We rarefied all libraries to 13,156 reads per sample for alpha and beta diversity analyses. We calculated alpha diversity metrics (Shannon index) using the 'estimate_richness' command in phyloseq. We visualised differences in Bray-Curtis distance among samples using non-metric multi-dimensional scaling (NMDS) ordination of bacterial communities using the R package vegan [33], setting k to 3 to yield a stress value of <0.17. We also used the vegan function 'vegdist' to calculate Bray-Curtis distance between bacterial communities to compare to genetic distances among samples, and the 'adonis' function to compare bacterial community structure among populations using permutational MANOVA. We used the R package gplots [34] to visualize the number of shared SVs among populations with a Venn diagram. To simplify the Venn diagram, we collapsed the two 'Tugela' sampling points into one group.

Analyses
Details of quality control and summary statistics carried out can be found in Supplementary Information. To assess levels of population differentiation, we calculated global F ST and pairwise F ST values between sites with corrections for null alleles using the ENA method of Chapuis and Estoup [37] with 1000 bootstrap replicates, and pairwise D [38] in Genodive v2.0b23. [39] We used ADZE v1.0 [40] to calculate site-level rarefied allelic richness (maximum g = 20) to study the effects of intraspecific genetic diversity at the population level on microbiome variation. We excluded Khubnam due to the low sample size (n = 9) compared with other sites (n = 20-28) to ensure sufficient power to detect differences among sites was retained. We created a pairwise matrix of Bruvo genetic distance [41] between all individuals in GenoDive. We also calculated individual-level heterozygosity (proportion of heterozygous loci in all loci successfully genotyped) using Genhet [42] in R to study its effect on associated microbial communities, as individual heterozygosity can be correlated with fitness.
We conducted principal coordinates analysis (PCoA) in GenAlEx v6.502 [49,50] using pairwise D values and the standardised covariance method. We ran an analysis of molecular variance (AMOVA) in Genodive using the Infinite Allele Model. We tested for the presence of isolation by distance among all locations and within population clusters. We conducted a Mantel test in ade4 in R [51] in R with 999 permutations to test associations between matrices of pairwise linearised F ST values (F ST /1-F ST ) and the logarithm of geographic distances. We used GENECLASS2 v2.0 [52] to identify putative first-generation migrants among sites [53,54].
We used INEST v2.1 [55] to identify the presence of two genetic signatures of bottleneck events; heterozygosity excesses in respect to allelic richness [56] and deficiency in M-Ratio values (mean ratio of the number of alleles to the range in allele size; [57]).

Individual level data
To examine the factors affecting Shannon diversity, we fitted a Gaussian Mixed effects model in the R package lme4 [58] with Shannon index as the response, Gosner stage and heterozygosity (specified as proportion of heterozygous loci, [59]) as predictors, and population ID as a random intercept (n = 126 data points for which there were both microbial (Shannon) and host genetic (proportion of heterozygous loci) data). Gosner stage and proportion of heterozygous loci were z-transformed prior to model fitting. The most complex model contained a proportion of heterozygous loci × Gosner stage interaction. We assessed significance of terms using likelihood ratio tests between models estimated using maximum likelihood. To quantify differences in mean richness among populations, we extracted the posterior modes of the population random effects from a Bayesian version of the best-supported model, specified in the R package MCMCglmm [60].
To examine the factors predicting microbial community structure (beta diversity), we calculated the withinpopulation 'divergence' metric of beta diversity using the microbiome package (Lahti et al. 2017) to be used as a response in a Gaussian mixed effects model. The community structure model contained Gosner stage, proportion of heterozygous loci and their interaction as fixed effects, and population ID as a random effect. We calculated the r 2 values of minimal models using the [61] method for mixed effects models, as implemented in the R package MuMIn [62]. We used a partial Mantel test to examine the correlation between genetic distances among A. hymenopus individuals (calculated from microsatellite data) and Bray-Curtis dissimilarities among their microbial communities (calculated from 16S amplicon data), while controlling for the effect of geographic distance. We used the 'mantel.partial' function in the vegan package, specifying the Pearson correlation statistic.

Population level data
We fitted population-level regressions (n = 6 sampling sites) to test the influence of allelic richness and Bd infection prevalence on alpha diversity. We did not fit both Bd prevalence and allelic richness in the same model to maximize residual degrees of freedom in the model. Significance of terms was assessed by likelihood ratio test.

Bd infection
The average Bd infection intensity across all individuals was 0.19 (±0.02) zoospores, with an average prevalence of 10.91 (±3.92)% across populations. There was relatively low infection prevalence and intensity at all sites, with highest infection rates at Tugela 2 and a complete absence of infection at Ribbon Falls (Table 2).
We observed significant differences in Shannon diversity index among sampling sites (Fig. 3). Tukelahed had the highest alpha diversity, whilst the values for the remaining five sampling locations were significantly lower (95% credible intervals of difference all >0; Fig. 3). There was no significant difference in the alpha diversity of the two Tugela sampling locations (95% credible interval of difference −0.92 to 0.05). NMDS ordination revealed distinct separation of bacterial community centroids corresponding to sampling site (Fig. 4). PERMANOVA of community distances supported this pattern, with significant separation among sampling sites (F 5,127 = 18.64, r 2 = 42.32%; p < 0.001). Shannon diversity of tadpole mouthparts was not affected by either Gosner stage of tadpoles (χ 2 = 2.61, p = 0.1), proportion of heterozygous loci (χ 2 = 2.19, p = 0.14), or their interaction (χ 2 = 0.19, p = 0.67).
Gosner stage significantly influenced within-population divergence in beta diversity (χ 2 = 5.01, p = 0.025; Fig. 5). There was no evidence of an effect of an interaction between Gosner stage and proportion of heterozygous loci on beta diversity, or proportion of heterozygous loci as a main effect (all p > 0.47). The best model explaining differences in beta diversity contained Gosner stage as the only fixed effect and had a marginal r 2 of 4%.
There was a positive correlation between host genetic distance calculated from microsatellites and microbial community distance when controlling for geographic distance (Partial Mantel Test r = 0.145, p < 0.001). There was no evidence that the Shannon diversity at a sampling site was affected by site-level allelic richness (F 1,4 = 0.95, p = 0.38) or Bd infection prevalence (F 1,4 = 1.17, p = 0.34).

Population genetics
Strong genetic structure was found in A. hymenopus, with F ST levels reaching 0.292 between Ribbon Falls and Nampolice (Table S2 and Supplementary Information). Ribbon Falls was highly differentiated from all sites, and in general sites from different rivers showed moderate levels of genetic differentiation, while sites within the same river showed low differentiation. STRUCTURE showed the optimum K to be two (Fig. 6), with Ribbon Falls forming a distinct population cluster and the remaining sites another. At K = 3, Vemvane and Nampolice also formed a separate cluster (Fig. 6). Hierarchical clustering analysis conducted on the population cluster of the K = 2 analysis (i.e., repeating the STRUCTURE run without Ribbon Falls) supported the presence of a Nampolice and Vemvane cluster. The PCoA and AMOVAs showed results consistent   with the STRUCTURE analysis (Fig. 7, Table S4), and together indicate that Ribbon Falls is a distinct population cluster, and Vemvane and Nampolice form a subpopulation cluster within the cluster of remaining sites.
Bottlenecks were detected at three sites: Heterozygosity excess in respect to allelic richness was found at Nampolice (p = 0.003), and deficiencies in M-ratios were found at Nampolice (p = 0.003), Vemvane (p < 0.001) and Ribbon Falls (p = 0.042).

Discussion
Population structure in A. hymenopus was strongly influenced by river network and geographic distance. However, at the landscape level, the microbiome of A. hymenopus is not structured according to host population genetic structure or river system. For example, Tukelahed did not have a similar microbial community structure to downstream sites Tugela and Tugela 2 (Fig. 4). Despite a lack of similarity between host population genetic structure and site-level microbiome structure, genetic distance among hosts is still significantly correlated with associated microbial community dissimilarity when controlling for geographic distance. This has not been previously shown in amphibians, but host genotype influences the gut microbiota of chickens [16], mice [63] and fish [7], and inbreeding and relatedness (although not population genetic structure) influence gut microbiomes of gopher tortoises (Gopherus polyphemus; [8]). Although not tested here, these patterns may be due to genetic variation in microsatellites reflecting variation in phenotypic traits that affect microbiome composition.  K h u b n a m N a m p o li c e V e m v a n e T u g e la 2 T u k e la h e d T u g e la B il a n ji l R ib b o n F a ll s K = 2 K = 3 Fig. 6 Bayesian population assignment analysis STRUCTURE plots showing membership coefficients (Q) for K = 2 and K = 3 genetic clusters in Amietia hymenopus. Vertical bars represent individuals, and colours represent a genetic cluster; the height of the bar represents the proportion of the genotype assigned to that cluster Microsatellites can be linked to genes coding for functional traits [64][65][66], including immune response genes in a number of species [67][68][69]. Variation in these genes could influence the composition of microbiomes by controlling the horizontal transmission of microbes from the environment. In A. hymenopus, the microbiome of the tadpole mouthparts is most likely transmitted from the environment as they are feeding or behaviourally selecting areas in the river. Behaviours can also influence symbiont transmission [70], and therefore there may be a link between genetic variation, behavioural variation, and microbe acquisition. That the influence of host genotype does not extend to landscape level microbiome structure may indicate that sitelevel environmental differences have more influence in structuring tadpole microbiomes, and host genotype simply exerts some selection pressure or behavioural influence on horizontal transmission. We detected a link between Gosner stage and mouthpart microbiome composition. Studies to date have predominantly used skin swabs to determine host microbiome of amphibians, as this is the region that becomes infected by Bd on postmetamorphic individuals. However, Bd primarily infects the mouthparts of larvae as keratin is localised to this area until the onset of metamorphosis when skin is remodelled to contain keratin. Although Bd infection intensities tend to increase as larvae develop due to increasing size [71], we show that alpha diversity (i.e., number of SVs) of the larval mouthpart microbiome does not increase as Gosner stage increases. This indicates that ontogenetic variation in microbiome is not simply an artefact of increasing mouthpart size, but rather it is the specific composition of the community that is changing as tadpoles develop. These microbiome changes may be related to phenotypic changes in mouthparts over development [24,72], which occur in other montane river species in the region during development as an adaptation to fast-flowing water [73]. Ontogenetic variation may also arise from changes in behaviour throughout larval development that means individuals are exposed to different environmental conditions and thus, microbial influences; however, this has not been studied in this species.
Despite these genotypic and ontogenetic influences on microbiome composition, populations shared a considerable core microbiome (Fig. 2), which has also been shown for other amphibian species [6,74], including among genetically isolated populations [75]. Our population genetics analyses detected low numbers of first-generation migrants between populations, with high-genetic differentiation given the relatively small geographic area studied. Together, these results indicate very low levels of dispersal among rivers, and in some instances, low dispersal among tributaries in the same river. This implies that core microbiome observed across the metapopulation of A. hymenopus has not arisen from contact between populations. A number of studies have now indicated that amphibians are able to select their microbiome based on functional traits [76,77]. This may be related to conserved host traits, such as anti-microbial peptide composition, or adaptive genetic loci. In addition, the river system and associated vegetation is remarkably homogenous on the uninhabited and unspoiled ecosystem of the Mont-aux-Sources plateau, with very low diversity of plants [78]. This may also in part explain the shared microbiome exhibited across populations. That said, there is likely to be fine-scale variation between environmental conditions in the rivers, which may be driving some variation in microbiome composition [6,79,80]. Microbial communities have been shown to vary between rivers even when colonisation substrate is controlled for [81]. Individuals from Nampolice and Tukelahed had the most distinct microbiome (Fig. 4), and the highest alpha diversity (Fig. 3), and this may be related to proximity from the source of the river. For example, Górniak et al [82] found bacterial diversity generally decreased as they sampled downstream from the source (glacial lake), and linked this to changes in environmental conditions and nutrient availability.
Historically there is wide annual and geographical variation in Bd infection prevalence in the Mont-aux-Sources region, ranging from~0 to 100% [23]. Although, we found uncharacteristically low infection loads in this study, our data correspond to historical patterns for the three sites that have undergone long-term monitoring. Ribbon Falls has historically low infection prevalence, followed by Bilanjil, and then Tugela, which shows the highest infection prevalence [23]. This may be driven by a range of environmental factors similar to those that drive microbiome diversity, including water flow, temperature, pH and depth [71]. Particularly low infection intensities such as those found here are, however, difficult to assess accurately using quantitative PCR. Given these low infection intensities, it is not possible to identify meaningful relationships between Bd infection and microbiome diversity, although this has been demonstrated in other studies ( [2,12,27,76,83,84]).
Ribbon Falls, Nampolice and Vemvane showed evidence of recent population bottlenecks, and these sites also have some of the lowest allelic richness. Considerable mortalities of A. hymenopus were seen at Ribbon Falls and Vemvane in September 2006 during routine monitoring (every 2-4 months) of Vemvane, Tugela, Bilanjil and Ribbon Falls ( [23]; C. Weldon, pers. obvs.). Diagnostic analyses identified Bd infections in 100% of the dead metamorphic A. hymenopus (n = 10) collected at each site. Further dead post-metamorphic individuals (n = 1-12) were encountered on most of the subsequent surveys, all infected with Bd [23]. Low genetic diversity in mortality-affected sites is unsurprising given the low level of migration from other sites. This presents a potential danger in bottlenecked locations struggling to recover evolutionary potential, impacting their ability to respond to stressors such as climate change, disease or local climatic events. Local chytridinduced mass mortality events may eliminate genetic diversity from the species as a whole due to the distinct population genetic structure and low migration. Interestingly, Bd infection prevalence at Ribbon Falls is generally low compared to other rivers in the region (this study and [23]), and Ribbon Falls has the lowest genetic diversity, which may be because of selection of Bd-resistant individuals. Addis et al. [85] showed montane boreal toads (Bufo boreas) with higher heterozygosity had higher Bd infection intensity, which they attribute to increased dispersal of heterozygous individuals leading to increased spread of Bd. Therefore, it is possible that the lack of immigrants to Ribbon Falls is also limiting the spread and maintenance of Bd at this site.

Conclusion
Here we have shown that the microbiome composition of amphibians is determined by genetic and ontogenetic variation, as well as geographic site, indicating within and among population predictors of microbiome diversity. Bd prevalence was low in the year we sampled, but we found lower genetic diversity and signals of genetic bottlenecks in sites where mass mortalities have been recorded in previous years, demonstrating the vulnerability to and persisting effects of disease in amphibian populations. Together these results demonstrate the key role of host genotype on microbial colonisation, which may be important in pathogen susceptibility in the amphibianchytridiomycosis system. Identifying individuals based on this characteristic may aid amphibian conservation efforts against this lethal pathogen.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons. org/licenses/by/4.0/.


Transforming Carcinogenesis is Essential Substrate Accumulation/Interactivity Dimension in Ubiquitin-Proteasomal Dysfunction

The constitutive susceptibilities to malignant transformation are well illustrated by the evolving dimensions of a dysfunctional ubiquitin proteasomal participation in generative reformulation of carcinogenetic pathways of interactivity. Further considerations indicate the transformation as derivative dysfunctional ubiquitin proteasomal system for further transformation. The multilayered complexities of such changes attest for redistribution of protein intermediates within various cellular subcompartments as signified by endosomal and nuclear systems of cooperative interchange. Dynamics formulation is itself a parametric measurement in substitution modelling of substrates as indeed projected within system profile for further change.Theincremental disjunction of such formulasareproposed alternative systems in further parahomeostatic dysregulation in carcinogenesis as a primary generative induction as signified by further developmental adaptation of such systems as DNA repair pathways.

Introduction
The ubiquitin proteasomal system is essential requisite schematic machinery in the operative dimensional functioning of whole pathway operability in homeostasis.Such system is hence potential dysfunctioning of the various mechanisms ranging from receptor turnover to endosomes, cell cycle dynamics, DNA repair processes. Ubiquilin4, a member of the UbL-UBA protein family, serves as an adaptor in the degradation of specific substrates via the proteasomal pathway and exerts anti-tumor potential in gastric cancer by p53dependent and p53independent regulation of p21 [1]. Ubiquitinspecific protease 15 is essential asanhuman papillomavirus16 E6interacting protein [2]. The ubiquitin and ubiquitinlike intermediates are hence a core system that operates also in induced dysfunction of cell disease and cancer.USP22 deficiency leads to myeloid leukemogenesis upon activation of oncogenic Kras activation through a PU1 dependent mechanism [3]. It is relative to such considerations that intense investigation has been devoted to modulate and suppress the dysfunctional proteasomal pathways to suppress carcinogenesis.

Dynamics of Operability
Dynamics of operability of the ubiquitin proteasomal system exhibit exquisite specificity for protein molecular targeting and hence are ideally constituted in the face of a broad range of protein substrates through posttranslational modification of numerous protein moieties. In such manner, evolving pathways can potentially incriminate the endless formation of monoubiquitylation conjugates and also polymeric chains of ubiquitin in the formulation of dysfunctional regulation of cellular pathways that otherwise are of significance in homeostatic control. HERC proteins exert ubiquitin ligase activity and are widely involved in various cancers [4]. Increments of operability include multilayered complexity in the system profile dysfunctions in cancer as wellillustrated by the various developmental pathways that operate in progression of carcinogenesis. Ubiquilin4 regulates the degradation of many proteins implicated and may induce cell cycle arrest and apoptosis via ERK activation and upregulation of cyclin D1 [5]. Carcinogenesis is a parahomeostatic series of mechanistic pathways that operates in terms of adaptor systems and as dysfunction adaptation to dysregulatory lesions within cellular pathways. Deubiquitinylase USP47 promotes RelA phosphorylation and survival ingastric cancer cells [6]. The incremental dimensions for polyubiquitylation include the creation also of polymer homotypic and also heterotypic ubiquitin chains that compete with ubiquitin and ubiquitinlikesystems for Glycineinduced attachments to the Lysine moiety of substrate proteins. SUMOylation is covalently modified on lysine residues of a target protein and regulates transcription, cell cycle, DNA repair and innate immunity [7].

Complexity
Many layers of complexity characterize single and multiple rounds of ubiquitylation in a manner that allows for crosstalk to ubiquitin like modulators such as NEDDylation and SUMOlylation of various protein substrates.Cullin Ring Ligases in particular are implicated in such crosstalk in a manner that implicates promiscuous interactivities of whole protein substrate pathways leading to carcinogenesis. STXBP4 regulates anaphase promoting complex or cyclostomes mediated p63 turnover and drives squamous cell carcinogenesis [8]. The development of parahomeostatic dysregulation is hence a prerequisite in the formulation of dimensional adaptation to cellular dysfunction. E3 ubiquitin ligase RAD18 and its downstream effectors of the Yfamily DNA polymerases participate in different phases of multistep carcinogenesis [9]. Developmental disorders implicate and further enhance the ubiquitin-proteasomal system within the various dimensions for complexity disorders. Histone chaperone ASF1A is physically associated with USP52 which in turn acts as a deubiquitinase in chromatin assembly with a potential role in breast carcinogenesis [10]. The ubiquitinspecific protease USP28 deficiency promotes breast and liver carcinogenesis as well as tumor angiogenesis in a HIF-independent manner [11]. The performance of incremental dysfunctionality is proof of principle towards the changes in profile interactivity of protein substrates that may compete with ubiquitin in the stabilization of suppressors and oncogenes such as p53 and HDM2 formulas. SUMO catalytic cycle includes maturation, activation, conjugation, ligation and demodification and is widely involved in carcinogenesis, DNA damage response, cancer cell proliferation, metastasis and apoptosis [12].

Derivative Phenomena
Derivativescope for multilayered complexities in ubiquitylation is an attempt for maladjustment in formulating the further comparative dimensions as projected by protein substrate interactions to active sites on other protein systems within the cell. The performance co-operativity of such system pathways hence augment the range of effects of operative ubiquitylation towards such profile systems as demonstrated by multi-compartment dysfunction within cells that subsequently undergo malignant transformation. The cullin family of proteins is implicated in the ubiquitin or mediated degradation of proteins, modulating cell proliferation, cellcycle control, migration, invasion and metastasis in tumor progression [13]. Wholegenome sequencing of colorectal cancer identified amplification of Ring finger protein 6, a RINGdomain E3 ubiquitin ligase, and plays a pivotal oncogenic role in colorectal tumourigenesis [14].
System adaptation in the formulation of dysfunctional profiles in ubiquitylation and of ubiquitylation-like formulas allow for extensive impact dynamics in the non-resolution of pathway generation and progression. Processivity systems permit the modulation of injury or lesion demarcation that comes to dominate in extensive manner the protein milieu of pathway progression in parallel modetoaffect multiple parahomeostatic systems of effect and transformation. Ubiquitination and ubiquitinlike post translational modifications regulate activity and stability of oncoproteins and tumor suppressors; USP18 is the major DUB that removes IFNstimulated gene 15 from conjugated proteins [15].

Maladaptation
Intense investigations have been carried out in profile formulation and require the distinctive global dimensions of cellular dysfunction that promotes accumulation of various projected moieties in manner cooperativity when the ubiquitinproteasomal system is inhibited.Such formulas of adaptive dysfunction include gene suppressors in a manner that is modeled by stabilization of the p53. Performance increments detail the constitutive and parahomeostasis that redefine the projected dynamics of turnover of proteins in terms that signify further operativetransformation of biochemical and biophysicaldimensions of maladaptation as well exemplified by nuclear repair attempts of damage.

Para Homeostasis
Cooperative parahomeostasis is hence complexitydriven in terms of the initiation and further constitutive progression of DNA damage and as abnormal cell compartmentalization of protein moieties.The distinctive compartmentalizationofsuch protein molecules incrementallyallow for permissive reconstitution as projected within schemes for further complexity formulation.The active participation of DNA damage and the accompanying systems in DNA repair systems are further dimensionalized as well developed maladaptive system characterization and recharacterization. RINGFinger protein 6 amplification activates JAK or STAT3 pathway by modifying SHP1 ubiquitylation and associates with poor outcome in colorectal cancer [16]. The pathways that formulate DNA damage are themselves modified in terms of milieu transformation within carcinogenetic initiation and progression.Ubiquitylation pathways initiated specifically by UBA6 set a suppressive barrier against critical steps of mammary carcinogenesis such as loss of polarity, anoikis resistance and epithelial mesenchymal transition [17]. Such pathway systems perform idealistic recapitulation within the broad range for competitive protein interactivity as induced by altered ubiquitin-proteasomal function and dysfunction.

Concluding Remarks
Derivative characterizations of inhibition and transformation of interactivities of the ubiquitin or proteasome is ideal model participation in the realization of system reformulation in a manner that promotes a series of adaptive transformations consequent to accumulative dynamics with substrateenhanced inhibition of the ubiquitylation of a broad range of protein metabolism and pathway crosstalk interactivity. MARCH8 belongs to a family of membrane associated RING-CH(MARCH) ubiquitin ligases and is associated with poor prognosis in non-small cell lung cancer patients [18]. Formulas are created and these generate various maladaptive transforming influences in the system profile projection of such global participation as carcinogenesis. Projected complexity systems well demonstrate the incremental formulas as pathway cross-talk and cooperative interaction as system modulatory pathways that conceive of generative transformation to carcinogenetic pathways for further transformation. Performance utilization pathways

3/5
induce enhanced process formulation that is mirrorimaged in the DNA damage to repair systems of adaptation and readaptation. The performance attributes for further transformation in terms of maladaptive change allow for incremental progression in cell cycling events and as selective participation that project with the onset impact of transforming parahomeostasis. Distributed compartmentalization of intermediate participants allows forcomplex microenvironmental modulation within transformed dynamics for carcinogenesis that is both adaptive and generative within the global cellular dimensions for projected complex interactivity.


text:
Microbiome Analysis Reveals the Presence of Bartonella spp. and Acinetobacter spp. in Deer Keds (Lipoptena cervi)

The deer ked (Lipoptena cervi) is distributed in Europe, North America, and Siberia and mainly infests cervids as roe deer, fallow deer, and moose. From a one health perspective, deer keds occasionally bite other animals or humans and are a potential vector for Bartonella schoenbuchensis. This bacterium belongs to a lineage of ruminant-associated Bartonella spp. and is suspected to cause dermatitis and febrile diseases in humans. In this study, we analyzed the microbiome from 130 deer keds collected from roe deer, fallow deer and humans in the federal states of Hesse, Baden-Wuerttemberg, and Brandenburg, Germany. Endosymbiontic Arsenophonus spp. and Bartonella spp. represented the biggest portion (~90%) of the microbiome. Most Bartonella spp. (n = 93) were confirmed to represent B. schoenbuchensis. In deer keds collected from humans, no Bartonella spp. were detected. Furthermore, Acinetobacter spp. were present in four samples, one of those was confirmed to represent A. baumannii. These data suggest that deer keds harbor only a very narrow spectrum of bacteria which are potentially pathogenic for animals of humans.

INTRODUCTION
Blood-sucking arthropods are vectors for numerous infectious agents in humans and animals and therefore of high interest in one health approaches. Deer keds (Lipoptena cervi) belong to the family of the louse flies (Hippoboscide) which are found in Europe, North America, and Siberia (Lindener, 1964). From August to November, the winged imagines (adults) fly to suitable hosts (mainly cervids) and lose their wings before they start to suck blood. They give birth to living larvae which fall to the ground as pupae and remain there until new imagines hatch to find new hosts (Haarløv and Haarlov, 1964). Figure 1 shows a scheme of the life cycle of deer keds.
It is unclear if the bites cause harm to the infested cervids: moose, which usually are highly infested with deer keds (Madslien et al., 2012), do not show worse indices of health compared to moose which live in deer ked-free areas (Paakkonen et al., 2012). However, hair loss in moose occurs when they are heavily infested with deer keds (Madslien et al., 2011). FIGURE 1 | Life cycle of deer keds. The winged adults fly to suitable hosts (mainly cervids) and lose their wings before sucking blood. They give birth to living larvae which fall to the ground as pupae and remain there until new adult deer keds hatch to find new hosts.
Bartonella spp. are Gram-negative, facultative intracellular bacteria which are usually transmitted by blood sucking arthropods and which can cause intraerythrocytic infections in their reservoir hosts (Dehio, 2005;Mändle et al., 2005;Maggi et al., 2011). The most frequently detected pathogen among in humans is B. henselae, the agent of cat-scratch disease. Cats are the reservoir hosts for these bacteria and pathogens are transmitted to humans by scratches or bites. Dogs can also become infected and develop endocarditis, fever of unknown origin and peliosis hepatis (Kitchell et al., 2000;Fenimore et al., 2011;Maggi et al., 2011;Drut et al., 2014). At least 37 Bartonella spp. are described which can infect a broad variety of mammals . B. schoenbuchensis was first isolated from the blood of wild roe deer in 1999 (Dehio et al., 2001). It belongs to a lineage of ruminant-associated Bartonella spp. comprising of B. schoenbuchensis, B. capreoli, B. chomelii, B. bovis, and B. melophagi (Engel et al., 2011). Reports of diseases associated with these bacteria are rare for animals and humans as chronic asymptomatic infections with long lasting bacteremia Abbreviations: A, Acinetobacter; B, Bartonella; L, Lipoptena; OUT, Operational Taxonomic Unit. are common for Bartonella spp. in their respective reservoir hosts (Rolain et al., 2001;Dehio et al., 2004). There is also evidence that these Bartonella spp. might cause endocarditis, fatigue, muscle pain and fever in their reservoir hosts or even in humans (Maillard et al., 2007;Maggi et al., 2009;Vayssier-Taussat et al., 2016). In this study we report a comprehensive microbiome analysis using next generation sequencing (NGS) to address further pathogenic agents in deer keds from Germany.

Sample Collection
The samples were collected between May and December 2017 from hunted roe and fallow deer at several locations in the federal state of Hesse, Germany, and nearby Karlsruhe in the federal state of Baden-Wuerttemberg, Germany. Deer keds were also sent from Ettlingen, Baden-Wuerttemberg, Germany and Wittstock, Brandenburg, Germany. Keds were collected in sterile, DNAfree tubes (Eppendorf, Hamburg, Germany) containing 70% ethanol. Whenever possible, muscle samples or full blood was collected from the hunted animals. No experimental procedures on animals or humans were performed. Only deer keds from hunted animals and from humans who sent them to us for diagnostic reasons were screened. For these procedures, there is no need for a permission from an ethics committee in Germany.

DNA-Extraction From Deer Keds, Muscle, and Full Blood
DNA extraction from deer keds and muscle samples of hunted deer was conducted as previously described . All deer keds were individually removed from the tubes with sterile forceps. Each deer ked was treated individually to prevent cross-contamination rinsed once in ethanol and twice in sterile water. DNA was extracted with the QIAamp DNA Mini Kit (Qiagen, Hilden, Germany) according to manufacturer's instructions. Grinding of deer keds and muscle samples was conducted with disposable, sterile mortars and pestles. DNA from full-blood was extracted using the DNeasy Blood and Tissue Kit (Qiagen) according to manufacturer's instructions.
The laboratories of the Institute for Medical Microbiology and Infection Control at the University Hospital of the Goethe University in Frankfurt (Germany) undergo a strict quality control management (DIN ISO 15189:2014 certificate, valid through January 2021). There was no increase of Bartonella or Acinetobacter positive cases during this study; therefore, the possibility of DNA contamination from non-study material is highly unlikely.

Microbiome Analysis of Deer Keds Using Next Generation Sequencing by Illumina Technology
The 16S rRNA gene of each deer ked DNA sample was amplified with primers for the V4 region (Caporaso et al., 2011) and analyzed as previously described (Regier et al., 2018). In brief, the V4 region amplification was done using Platinum SuperFi PCR Master Mix (Thermo Fisher Scientifc, Carlsbad, U.S.A.). Thermocycler conditions were used as described previously (Regier et al., 2018). PCR products were purified using AMPure XP DNA beads (Beckman Coulter, Brea, U.S.A.) before running the index and adapter ligation PCR with the Nextera XT Index Kit v2 Set A and B (Illumina, San Diego, U.S.A.) as described by the manufacturer. Quality controls of libraries were done using the Qubit 2.0 instrument (Thermo Fisher Scientific) and the 2100 Bioanalyzer instrument (Agilent Technologies, Santa Clara, U.S.A). Samples were diluted, pooled, spiked with an internal control (15% PhiX) and paired-end sequenced on the MiSeq Illumina platform using a flow cell with V2 chemistry (500 cycles). Negative controls were performed using pure water and elution buffer. In addition, microbial mock communities (Zymo Research, Irvine, California, USA) were run along as a standard and as a quality control for determining contamination bias from DNA extraction.

Bioinformatic Microbiome Analysis Workflow
The bioinformatics analysis was performed as previously mentioned (Regier et al., 2018). Briefly, the paired end reads were joined and the primer sequences were removed. Reads with ambiguous base calls or with homopolymers longer than eight nucleotides were removed and duplicate sequences were merged and aligned against the SILVA-bases bacterial reference alignment (Quast et al., 2013). Using the Mothur implementation of the uchime algorithm, chimeric reads were removed, taxonomy was assigned and non-bacterial reads were removed. OTUs were created using Mothur and the taxonomy was reassigned to the ladder. In preparation for the analysis with Qiime, a phylogenetic tree and an OTU table in biom format was created. Subsequently, the alpha-diversity analysis and the taxa summary plots were created using the Qiime core diversity analysis script.
Confirmation and Species Determination of Bartonella spp. by 16S-23S-rDNA-ITS-PCR and rpoB-PCR Bartonella positive deer keds and corresponding muscle and full-blood samples, if available, were screened for Bartonella DNA. To detect Bartonella spp. specific DNA in deer keds, two PCRs were conducted: a 16S-23S-rDNA Internal Transcribed Spacer (ITS)-region-PCR (Cherry et al., 2009) and a PCR detecting the sequence for the rpoB gene (encoding the βsubunit of the bacterial RNA polymerase) was performed as previously described (Oksi et al., 2013). Both PCRs were conducted with the Platinum Taq Polymerase-Kit (Invitrogen, Schwerte, Germany). All PCR primers are listed in Table 1. Positive (B. henselae Houston ATCC 49882) and a negative (water) controls were always included. DNA was amplified in a Biometra T3000 thermocycler. Products were separated on an agarose gel, ethidium bromide stained and visualized under UV light. All PCR-products were sequenced (GATC, Konstanz, Germany). Sequences were analyzed using Chromas software (Technelysium, Version 2.6, South Brisbane, Australia) and compared to Bartonella spp. strains deposited in the NCBI databank using BLAST online tool to distinguish them on the species level.

Nucleotide Sequence and Phylogenetic Analyses
Type alleles of the deer ked-derived partial rpoB fragments (Bs_GER_A GenBank: MH598359, Bs_GER_B GenBank: MH598360, and Bs_GER_C GenBank: MH598361) were used to search for similar sequences in nucleotide databases by Standard Nucleotide BLAST (BLASTN) at NCBI (https://www.ncbi.nlm. nih.gov) in June 2018. Search criteria were formatted with Expect min as 96 % and Expect max as 100 % to exclusively focus on close by entries within the ruminant lineage Bartonella species. Sorting was based on query coverage. BLASTN hits for all three rpoB type alleles were pooled into a single dataset, including removal of duplicate and triplicate entries. Phylogenetic analyses were performed using Molecular Evolutionary Genetics Analysis (MEGA) 6.06-mac (www.megasoftware.net/). To this end, the sequences were first aligned with ClustalW. Next, the partial rpoB fragments were either trimmed to a 406 bp fragment corresponding to a B. henselae Houston-1 rpoB fragment (AF171070) between nucleotide positions 1710 (TCGT...) and 2115 (. . . TCCA) or to a 285 bp fragment corresponding to a B. henselae Houston-1 rpoB fragment (AF171070) between nucleotide positions 1747 (ATTG....) and 2031 (. . . AGTA). The first rpoB sequence trimming corresponds to the maximum length fragment which can be obtained with rpoB-specific PCR primers prAPT0244 and prAPT0245 ( Table 1). The second trimming corresponds to the maximum length fragment that is available for all the nucleotide database entries identified according to above BLASTN search criteria.
Confirmation and Species Determination of Acinetobacter spp. by Oxa51-PCR Acinetobacter spp. positive deer keds and corresponding muscle and full-blood samples, if available, were screened for Acinetobacter DNA. To detect A. baumannii-specific DNA, the gene for the A. baumannii-specific carbapenemase bla Oxa−51 was detected (Woodford et al., 2006), however, this chromosomally encoded carbapanemase does not contribute to carbapenem resistance of A. baumanni because it is not expressed. Primers are listed in Table 1. Positive (A. baumannii, patient isolate) and a negative (water) controls were always included.

Sample Collection
130 deer keds were collected from 39 roe deer (n = 109), 8 fallow deer (n = 13) and 2 humans (n = 8). While the deer keds collected from animals had already shed their wings and a blood meal, the samples taken from humans still had their wings and did not start to feed on their hosts. Whenever possible, blood, or muscle samples were collected from the host animals. Full blood of 5 roe deer and muscle samples of 2 roe deer were taken. Locations of ked collections are given in Figure 2A,C. The number of sampled animals is summarized in Figure 2B.

Next Generation Sequencing of Deer Keds for 16S rRNA Microbiome Analysis
To date, a broad and in detail microbiome investigation of whole deer keds has not been conducted. It has been reported previously that deer keds possibly also act as vectors for pathogenic bacteria (e.g., B. schoenbuchensis, (Vayssier-Taussat et al., 2016)). Hence, we were interested in identifying the microbial composition of deer keds sampled from different hosts throughout Germany. In total, 130 deer ked samples were paired-end sequenced on the MiSeq Illumina platform, each resulting in a minimum sequencing depth of 5,000 reads. Sequences of samples with <5000 reads were excluded.
The alpha diversity of deer keds sampled from humans and fallow deer reveals a higher number of OTUs compared to roe deer. This demonstrates higher species richness in deer keds sampled from humans and fallow deer (Figure 3).
To examine the microbial taxonomic distribution of sampled deer keds, cumulative bar charts comparing relative family abundances were created. Comparing all three groups, we observed two dominant OTUs which are present, Enterobacteriaceae and Bartonellaceae. Bartonellaceae is found in a higher abundance in the roe deer group (∼75%). The group of fallow deer (∼40%) and humans (∼20%) show a lower abundance of Bartonellaceae. The family of Enterobacteriacea, which was later identified as Arsenophonus lipopteni a known obligate intracellular symbiont of the deer ked (Nováková et al., 2016) was most dominantly abundant in the human group (∼70%) followed by fallow deer (∼45%) and roe deer (∼20%) ( Figure 4A). As seen in Figure 4B, which demonstrates the variation in relative abundances of the top 17 OTUs between the three groups, the group of fallow deer revealed a higher abundance of Ruminococcaceae, Veillonellaceae, and Prevotellaceae, which all are representatives of the oral or intestinal microbiome. Interestingly, Acinetbacter spp. was also in four deer ked samples, one which was identified as A. baumannii.

Deer Keds
Confirmatory PCRs on 95 deer ked samples revealed the presence of B. schoenbuchensis DNA in 93 samples. Three unique rpoBalleles (Bs_GER_A GenBank: MH598359, Bs_GER_B GenBank: MH598360, Bs_GER_C GenkBank: MH598361) were detected. The rpoB allele "Bs_GER_A" was found in one deer ked from Karlsruhe, Germany, harvested from a fallow deer (Dama dama). The highest BLASTN sequence identity score (query coverage sorting) of the Bs_GER_A allele was 99.2 % (377/380 bp) with the rpoB fragments obtained from cattle blood samples in Spain (KM215709) or from elk blood samples in USA (HM167505). Bs_GER_A allele also had high BLASTN sequence identity scores to short rpoB fragments obtained from moose blood samples in Finland, e.g., 100% to KU254139. Phylogenetically the Bs_GER_A allele clustered with rpoB sequences retrieved from ruminant blood samples in Finland, Japan, Spain, and USA ( Figure 5). The rpoB allele "Bs_GER_B" was detected in nine deer keds from Giessen, Germany and Frankfurt, Germany. These deer keds were either harvested from a fallow deer (Dama dama; n = 7) or from a roe deer (Capreolus capreolus; n = 1). The highest BLASTN sequence identity score (query coverage sorting) of the Bs_GER_B allele was 100 % (380/380 bp) with several rpoB fragments obtained from deer ked samples in Poland (e.g., MF580675). Phylogenetically, the Bs_GER_A allele also clustered strongly with short rpoB sequences retrieved from deer ked samples in Norway, e.g., JN990612 (Figure 5). The most common rpoB allele was "Bs_GER_C", which was found in 83 deer ked samples and from all sampling sites. All of these deer keds were harvested from a roe deer (Capreolus capreolus). The highest BLASTN sequence identity score (query coverage sorting) of the Bs_GER_C allele was 100 % (387/387 bp) with the rpoB fragments obtained from a roe deer (Capreolus capreolus) blood B. schoenbuchensis type strain R1 isolate in Germany (AY167409) or from a human blood B. schoenbuchensis strain MVT06 isolate in France (HG977196). These similarities were also reflected in the phylogenetic clustering of the Bs_GER_C allele (Figure 5). No Bartonella spp. were detected in deer keds collected from humans. An Acinetobacter spp. OTU (operational taxonomic unit) was found in four deer ked samples. In one sample, the presence of A. baumanii DNA was confirmed by detecting the A. baumannii specific OXA 51 gene. Table 2 shows a summary of the allele distribution among all samples as well as the hosts and sampling sites.

Deer Blood
Blood from five roe deer (Capreolus capreolus) was available. The corresponding keds were B. schoenbuchensis-positive. DNA extracted from the blood was analyzed for the presence of Bartonella DNA by rpoB-specific PCR and Sanger-sequencing. One blood sample was positive with 100 % sequence identity with the Bs_GER_C rpoB allele. One sample showed an rpoB allele (GenBank: MH598362), which was different from deer ked "Bs_GER_A, " "Bs_GER_B, " and "Bs_GER_C" rpoB alleles (max ID). The highest BLASTN sequence identity score (query coverage sorting) of the unique deer blood allele was 100% (406/406 bp) with the rpoB fragment of the B. capreoli type strain IBS193. Therefore, this blood sample was positive for B. capreoli, not for B. schoenbuchensis. Three blood samples remained negative. A. baumannii was not detected in any of the blood samples, analyzed with OXA 51 gene-specific PCR.

Deer Muscle
Muscle samples from two roe deer (Capreolus capreolus), whose keds were B. schoenbuchensis-positive, were also analyzed for the presence of Bartonella DNA by rpoB-specific PCR and Sangersequencing. One sample remained negative while the other sample was positive for B. schoenbuchensis DNA with 100 % identity to the most prevalent rpoB allele "Bs_GER_C" detected in the deer keds.

DISCUSSION
Blood-sucking arthropods are vectors for many human and animal pathogenic bacteria. Whereas a plethora of data regarding transmission of pathogens to animals and humans by ticks and fleas is available (Regier et al., 2016), virtually nothing is known about pathogen transmission by deer keds. Based on our recent experiences in performing tick metagenomics (Regier et al., 2018), where we identified in total six potentially pathogenic genera being present within the ticks, we expected a broad variety of human or animal pathogenic bacteria in these deer keds. Furthermore, Borrelia burgdorferi DNA and Anaplasma phagocytophilum DNA has already been detected in deer keds in earlier studies (Buss et al., 2016). However, we found a surprisingly low diversity of bacteria within these insects suggesting that deer keds are not transmitting a broad spectrum of pathogens. The deer ked microbiome mainly consisted of two OTUs: Arsenophonus spp. and Bartonella spp.. Arsenophonus lipopteni is a well-known endosymbiont in deer keds with unknown biological function (Nováková et al., 2016). B. schoenbuchensis DNA was previously found in moose, roe deer, red deer, and cattle so probably these ruminants represent the reservoir hosts for these species (Chang et al., 2000;Bermond et al., 2002;Rolain et al., 2003;Maillard et al., 2004;Adamska, 2008;Duodu et al., 2013;Welc-Faleciak et al., 2013;Korhonen et al., 2015). Lipoptena cervi is suspected to act as the main vector for B. schoenbuchensis, it was found by cultivation and via molecular methods in adult L. cervi (Dehio et al., 2004;Matsumoto et al., 2008;Duodu et al., 2013;Bruin et al., 2015;Korhonen et al., 2015;Szewczyk et al., 2017), L. mazamae (Reeves et al., 2006) and feeding ticks (Matsumoto et al., 2008). Furthermore, B. schoenbuchensis was shown to colonize the midgut of deer keds (Dehio et al., 2004). Replication in the arthropod is a crucial prerequisite for deer keds to act as a natural reservoir host or a natural vector. Several studies showed the presence of Bartonella DNA in deer ked pupae (Duodu et al., 2013), in pupae and adult winged deer keds (Korhonen et al., 2015) and of B. schoenbuchensis DNA in winged and wingless deer keds and in larvae (Bruin et al., 2015). Also, no Bartonella spp. were cultured from moose with no deer ked infestation (Duodu et al., 2013). All these findings make it very likely that L. cervi is the vector for B. schoenbuchensis and that the bacteria can be transmitted transstadially. One deer ked sample was positive for B. capreoli but the corresponding deer ked was not. This could lead to the conclusion, that B. capreoli is not transmitted by deer keds. Epidemiologically, our data suggest that B. schoenbuchensis alleles derived from fallow deer and roe differ as the rpoB allele Bs_GER_B is nearly exclusively seen in fallow deer and Bs_GER_C in roe deer. Reasons for this fact remain speculative but host specificity is common among Bartonella spp. (Dehio, 2005;Mändle et al., 2005;Maggi et al., 2011) and this finding might demonstrate the process of host adaption of B. schoenbuchensis by so far unknown molecular mechanisms. It might be speculated that fallow deer had to less time to adapt to B. schoenbuchensis Bs_GER_C and roe deer to B. schoenbuchensis Bs_GER_B, respectively. Although no Bartonella spp. were detected in deer keds collected from humans, there is a possibility for humans to be infested with these bacteria since Bartonella spp. were detected in a huge number of ked samples. To date, it is unclear if and which diseases can be caused by B. schoenbuchensis in humans but it is suspected to cause the so called "deer ked dermatitis." Persistent, therapy resistant, pruritic papules can form one to 24 h after deer ked contact and it was shown that immunologic mechanisms are probably involved in the pathogenesis (Rantanen et al., 1982). Dehio et al. (2004) proposed B. schoenbuchensis as a cause of the dermatitis because of the similarity to the primary manifestation of cat scratch disease caused by B. henselae but clear evidence for this is still missing. Other ruminant-associated Bartonella spp. have also been shown to cause diseases in humans and animals. The role of B. bovis as the causative agent of endocarditis in cows has already been reported twice (Maillard et al., 2007) but a B. bovis bacteremia had no effect on milk production or reproduction in cattle (Maillard et al., 2006). B. melophagi was isolated from the blood of two sick women suffering from fatique, muscle weakness, muscle pain and fever, but its role in the pathogenesis of these symptoms remains unclear . In another study, B. schoenbuchensis was isolated from the blood of a patient again suffering from fatigue, muscle pain and fever. The patient had a history of tick bites and was seronegative for Lyme borreliosis (Vayssier-Taussat et al., 2016). To analyze whether deer ked transmitted B. schoenbuchensis might contribute to the pathogenesis of deer ked dermatitis or to those unspecific entities attributed to ruminant-associated Bartonella spp., it would be necessary to perform, e.g., studies in which the presence of anti-B. schoenbuchensis-antibodies in deer ked-exposed patients would be systematically analyzed. However, serological tools to perform such surveys are not available and, moreover, no cut-off values for B. schoenbuchensis serology have been determined.
In 2017, the World Health Organization listed A. baumannii as one of the top pathogens for which new antibiotics are urgently needed (WHO., 2017). It is known that A. baumannii is present in livestock (e.g., chicken and geese) and in wild storks (Wilharm et al., 2017). Acinetobacter DNA was also found in ectoparasites of domestic animals (Kumsa et al., 2012). Moreover, A. baumannii DNA is present in up to 21% of body lice (La Scola and Raoult, 2004). Therefore, transmission of this emerging pathogen might be promoted by various arthropod species. Although we detected non-baumannii Acinetobacter spp. in only three specimens and A. baumannii only once in a deer ked, it can nevertheless be discussed that deer keds contribute to the distribution of A. baumannii in animals and humans. In conclusion, L. cervi should be considered as a highly likely vector for B. schoenbuchensis and a potential vector for Acinetobacter spp.. Because of the fact that symptoms attributed to B. schoenbuchensis infections are of limited scientific evidence, it might be worth in future to systematically analyze whether B. schoenbuchensis transmitted by deer keds contributes infectious diseases in animals and humans.

AUTHOR CONTRIBUTIONS
VK performed the experimental design. YR performed sample collection, DNA extraction and PCRs. KK, MW, and TH performed next generation sequencing and bioinformatic analysis. AP performed phylogenetic analysis. YR, KK, MW, AP, SG, TH, and VK performed writing of the manuscript. All authors read and approved the final manuscript.

FUNDING
This work was partially supported by a grant from the Bayer Animal Health Company (VK), Leverkusen, Germany and by the Robert Koch-Institute, Berlin, Germany (Bartonella consiliary laboratory, 1369-354, VK) and the Deutsche Forschungsgemeinschaft (DFG Forschergruppe 2251, Project B3 (VK, SG). The work of TH was supported by grants of the Deutsche Forschungsgemeinschaft SFB-TR84 project B01 (TRR 84/2 2014) and B08 (TRR84/3 2018). KK and TH were supported by BMBF CAPSyS project (01ZX16004C). The research in the laboratory of AP is financially supported by Academy of Finland (Academy Project, No 295296, 2016-2020, Sigrid Jusélius Foundation (Junior group leader grant, 2015-2018) and University of Turku, Finland. The deer ked microbiome project in the laboratory of AP was in part supported by a 2013 Jenny and Antti Wihuri Foundation grant.