Seasonal Changes in the Endosymbiotic Consortia of Aphids from the Genus Cinara

Buchnera aphidicola is the primary endosymbiont of aphids with which it maintains an obligate mutualistic symbiotic relationship. Insects also maintain facultative symbiotic relationships with secondary symbionts, and Serratia symbiotica is the most common in aphids. The presence of both symbionts in aphids of the subfamily Lachninae has been widely studied by our group. We examined two closely related aphids, Cinara tujafilina and C. cedri in the present study. Even though both B. aphidicola strains have similar genome sizes and gene contents, the genomes of the two S. symbiotica strains were markedly different. The SCc strain has the smallest genome known for this species, while SCt possesses a larger genome in an intermediate stage between the facultative S. symbiotica of Acyrthosiphon pisum (SAp) and the co-obligate S. symbiotica SCc. Aphids are vulnerable to high temperatures. Previous studies indicated that S. symbiotica SAp confers resistance to heat-shock stress. In order to clarify whether S. symbiotica strains from genus Cinara also play a role in heat stress protection, we performed a quantitative determination of the consortium Buchnera/Serratia from two geographically close populations, each of which belonged to the Cinara species examined, over two years in natural environments. We found no variation in the consortium from our C. cedri population, but a positive correlation between both endosymbiont densities and average daily temperatures in the C. tujafilina population. Even though S. symbiotica SCt may retain some protective role against heat stress, this does not appear to be due to the release of protective metabolites by cell lysis.

Mutualistic symbiosis between intracellular bacteria and insects is widespread in nature and appears to be one of the keys to the evolutionary success of this group of organisms (34). In most known cases, bacteria supply essential nutrients that are absent or scarce in unbalanced host diets (e.g., plant sap, vertebrate blood, and cereals). Based on their essentiality for host survival and fitness, two classes of insect endosymbionts have been identified, i.e., obligate or primary (P) and facultative or secondary (S).
Obligate mutualists are confined within specialized cells called bacteriocytes, and are transmitted strictly vertically from the mother to offspring. As expected from this lifestyle, they maintain highly constrained relationships with their hosts (51). In aphids, the P-endosymbiont Buchnera aphidicola has established long-term obligate mutualistic symbiosis (31). It provides essential amino acids in order to complement the phloem diet of its host, which is poor in nitrogen compounds (8). This mutualistic relationship has been supported by genomic (7,22,28,41,46,47,50), transcriptomic (16), and proteomic (42) analyses of several B. aphidicola strains belonging to aphids from different subfamilies.
Facultative endosymbionts are not essential for host survival, and beneficial, harmful, or neutral effects have been described (12,39,40). They have established a more recent association with the host than P-endosymbionts. Although they are also vertically transmitted, horizontal transfer events between hosts have also been reported (18,39,43). They may live in different locations inside the host, such as sheath cells, the hemolymph, secondary bacteriocytes, or surrounding primary bacteriocytes. Three main S-symbionts have been described in aphids, i.e., Hamiltonella defensa, Regiella insecticola, and Serratia symbiotica (32), with the latter being the most widely distributed among different aphid lineages (22,43,53). Most experimental studies have been performed on the pea aphid Acyrthosiphon pisum (subfamily Aphidinae), in which the three aforementioned S-symbionts have been found (44), and different beneficial effects have been proposed. All three common facultative symbionts of A. pisum are known to be involved in protection against parasitoid wasps (37,45), while R. insecticola provides resistance to the fungal pathogen Pandora neoaphidis (45). R. insecticola and H. defensa may also contribute to host plant specialization (14,18,49), while the main benefit provided by S. symbiotica is protection against heat stress (6,30,44).
The genome sequencing of S. symbiotica SAp (from A. pisum, subfamily Aphidinae) corroborated the facultative status of this symbiont (4). S. symbiotica is prominent in members of the subfamily Lachninae (3,22) and has been characterized in detail in two aphid species from the genus Cinara, the cedar aphid C. cedri (Mimeur, 1936) and thuja aphid C. tujafilina (del Guercio, 1909) (23,29). Although both Cinara species are clustered as sister clades, distantly related with the Aphidinae, phylogenetic analyses on S. symbiotica showed that SAp and SCt (from C. tujafilina) belong to the same clade (clade A), while S. symbiotica from other Lachninae (including SCc from C. cedri) belong to a different one (clade B; 22,29). In C. cedri, the genome sequencing of S. symbiotica SCc revealed that it has taken over some of the roles that B. aphidicola plays in other aphids. Thus, it is necessary to complement B. aphidicola BCc for the synthesis of tryptophan (15) and it is responsible for the synthesis of vitamins and cofactors (23). Therefore, both endosymbionts are considered to be co-primary. On the other hand, the genome sequencing of S. symbiotica SCt and its comparison with S. symbiotica SCc and the facultative S. symbiotica SAp revealed that it is in an intermediate stage between them (29). Moreover, while S. symbiotica SCt is rod-shaped, may be found extracellularly, and possesses a large genome (2.5 Mb), resembling the facultative strain S. symbiotica SAp (4), S. symbiotica SCc is always confined inside bacteriocytes, has a pleomorphic shape similar to B. aphidicola, and has the smallest genome known for this species, with a gene content similar to that of other P-endosymbionts (23). Nevertheless, the pathway for the biosynthesis of riboflavin was shown to be lost in the ancestor of B. aphidicola from the Cinara species, causing the fixation of S. symbiotica as an obligate endosymbiont (29). Differences between both S. symbiotica strains are marked, considering B. aphidicola strains BCc and BCt have similar genome sizes and gene numbers (24).
One important factor in determining the geographical distribution of a species is the range and variability of temperatures tolerated by it. Since they are poikilothermic, insects are very sensitive to extreme temperatures. Temperature not only has direct effects on insect hosts, but also on their symbionts, as demonstrated in field and laboratory experiments (9,13,33,44,48). In the case of aphids, in addition to their difficulty tolerating high temperatures (6,36), the thermal sensitivity of their obligate endosymbiont B. aphidicola further constrains the tolerance of the host (51). Several hypotheses have been proposed to explain the high susceptibility of the endosymbionts to heat stress, and are mainly related to the special genomic features of these microorganisms (51). They typically have a very high AT content, leading to more thermosensitive DNA molecules (52), but may also affect the stability of structural RNAs. The accumulation of non-synonymous substitutions in endosymbiont genomes due to their accelerated mutation rates may cause higher thermal susceptibility in the encoded proteins. Additionally, their reduced genomes have lost most genes encoding cell envelope proteins, potentially resulting in fragile cells.
Studies conducted on A. pisum demonstrated that heat stress reduced the number of bacteriocytes, causing important decreases in B. aphidicola densities (19,30,36), which has marked consequences on host fitness. However, the concomitant presence of S. symbiotica (and, to a lesser extent, H. defensa) may improve aphid fitness under heat stress (2,6,17,30,44). Previous studies suggested that the facultative symbiont enhances the preservation of bacteriocytes and, thus, the survival of B. aphidicola at high temperatures (30,44). Consequently, S. symbiotica infections are highly prevalent in pea aphid populations after periods of summer heat or in hot desert sites (30,51).
If the vulnerability of obligate endosymbionts to heat stress limits the tolerance of the host to temperature changes (51), the transition of S. symbiotica from a facultative to an obligate lifestyle may have reduced its ability to protect the aphid host from heat stress. The different level of mutualistic integration of the two S. symbiotica strains from Cinara aphids offers a good opportunity to test this hypothesis. In order to investigate whether S. symbiotica SCc and SCt play a role in heat stress protection, and if their presence affects B. aphidicola density, we performed a quantitative determination on the consortium Buchnera/Serratia from a single population of each Cinara species under study over two years in natural conditions in the Mediterranean area of Valencia (Spain).

Materials and Methods
Insect rearing and sampling C. cedri specimens were collected from a permanent population that was established from individuals collected in 2011 (23) and maintained since then on two cedar trees of the species Cedrus atlantica glauca in the facilities of the ICBiBE at the University of Valencia (Paterna, Valencia, Spain. 39° 30' 58" N, 0° 25' 20" W; 24). C. tujafilina aphids were collected on April 14, 2012 from thuja bushes of the species Platycladus orientalis from Paiporta (Valencia, Spain. 39° 25' 52" N, 0° 24' 59" W), and were used to infect two thuja bushes grown in 5-L pots (20×20×19 cm 3 ), which were initially maintained in Paiporta and transferred to the facilities of the ICBiBE on June 17, 2013. Trees and bushes were searched for adult aphids on the plant stem and leaves twice a week in 2012, and weekly in 2013. When adult aphids were scarce, we also included L4 instar insects in the sampling. Daily hourly temperatures in the closest meteorological station (8414A Valencia/Aeropuerto, 39° 29' 06" N, 0° 28' 29" W) were obtained from the Spanish National Meteorology Agency (AEMET; http://www.aemet.es). Since samples were always taken early in the morning, the temperatures considered for each sample were those recorded the day before. The temperature regime was characterized using the following daily parameters: maximum temperature (T max ), minimum temperature (T min ), average temperature (T ave ), temperature oscillation (T osc ), and the standard deviation of temperatures (SD T ).

Total DNA extraction and quantification
Total insect DNA ( T DNA) was extracted immediately after sampling using the method of Latorre et al. (26) with slight modifications. Approximately 6.5 mg of insects was gently homogenized in 160 μL buffer I (10 mM Tris-HCl [pH 7.8]; 60 mM NaCl; 5% sucrose; 10 mM EDTA) at 4°C. After the addition of 200 μL of buffer II (300 mM Tris-HCl [pH 8.0]; 1.25% SDS; 5% sucrose; 10 mM EDTA), the mixture was incubated at 65°C for 30 min, neutralized with 60 μL 3 M potassium acetate (pH 5.0), and kept at -20°C for 20 min. T DNA was concentrated by precipitation using standard protocols, resuspended in ultrapure water, and stored at -20°C until its use. The concentration and quality of T DNA were measured using a PicoGreen dsDNA Quantification Assay (Invitrogen [Thermo Fisher Scientific], Waltham, MA, USA).

Quantitative real-time PCR (qPCR)
The abundance of B. aphidicola and S. symbiotica in the two Cinara species under study was determined by qPCR on a LightCycler 2.0 with 20-μL capillaries using the LightCycler FastStart DNA Master PLUS SYBR Green I Kit (Roche Diagnostics) according to the manufacturer's instructions.
Multiple alignments were performed in ClustalW (25; http:// www.clustal.org) to detect adequate sequence regions for the design of specific primer pairs to the single-copy genes groEL from B. aphidicola and atpD from S. symbiotica (Table 1). The gene atpD was selected for the quantification of S. symbiotica because it is not present in any of the B. aphidicola strains under study (24,41). The sequences of the genes under study were retrieved from GenBank (Acc. No. for the complete genomes of B. aphidicola BCc, NC_008513.1; B. aphidicola BCt, CP001817.1 and S. symbiotica SCc, CP002295.1. Acc. No. of the scaffolds of the S. symbiotica SCt genome containing the genes atpD, FR904236.1, and groEL, FR904230.1). Primer pairs were designed with IDT's Oligo Analyzer 3.1 (http://eu.idtdna.com/analyzer/Applications/OligoAnalyzer). By taking advantage of the sequence similarity between the groEL orthologues from B. aphidicola BCc and BCt, a single primer pair was designed for both bacterial strains. The specificity of each primer set was checked by BLAST (1, http://blast.ncbi.nlm.nih.gov/ Blast.cgi) and empirically by gel electrophoresis, as recommended by the MIQE guidelines (5).
We used recombinant plasmid DNA as quantification calibrators. Each primer set was used for conventional PCR amplification with the KAPATaq DNA Polymerase Kit (Kapa Biosystems, Woburn, MA, USA) using insect T DNA as a template. The thermal cycling protocol was as follows: initial denaturation at 95°C for 2 min, followed by 30 cycles at 95°C for 30 s, at 60°C for 30 s, and at 72°C for 1 min; a final extension at 72°C for 2 min. Amplicons were purified with the High Pure PCR Product Purification Kit (Roche), and cloned into a pGEM-T Easy Vector (Promega, Madison, WI, USA). The recombinant plasmids pGEM-groELBCc, pGEM-groELBCt, pGEM-atpDSCc, and pGEM-atpDSCt were purified using the High Pure Plasmid Isolation Kit (Roche). Plasmid concentrations were measured using a PicoGreen dsDNA Quantification Assay (Invitrogen). Calibration curves were obtained according to Lee et al. (27). The regression lines for the standard curves had a mean squared error <0.1.
qPCR experiments were performed in a total volume of 10 μL, using 0.5 μM of each primer, and 1 μL of T DNA as the template.
Each experiment was performed simultaneously for both target genes, in two biological replicates, and in three technical replicates containing 2.0, 1.0, and 0.5 ng μL -1 of T DNA, respectively. Each technical replicate was analyzed twice. Following MIQE guidelines (5), qPCR reactions containing water or the corresponding recombinant plasmid instead of T DNA were also performed for each primer pair as negative and positive controls, respectively. qPCR reaction conditions were as follows: one cycle of 95°C for 10 min, followed by 35 cycles of 95°C for 5 s, 60°C for 6 s, and 72°C for 7 s. At the end of each run, the fidelity of the amplification was checked through a melting-curve analysis for each amplicon, with a temperature gradient of 0.1°C s -1 from 70 to 95°C. The technical replicates for each dilution had a mean squared error ≤ 0.2. The amplification results were examined by electrophoresis on 1.4% agarose gels. The efficiency of each biological replicate was calculated after each qPCR run, and only replicates with efficiencies equal to or greater than 1.7 were taken into consideration. Therefore, one C. cedri sample (June 10, 2013) and five C. tujafilina samples (April 29, July 29, October 31 in 2012; June 24 and July 1 in 2013) were excluded from subsequent analyses (black dots in Fig. 1).

Statistical Analysis
In order to rule out that a putative effect of temperature on the variable examined was spurious and actually due to aphid biomass, we initially examined the relationship between aphid weight (mg) and each of the temperature parameters under study. We found no correlation between the aphid biomass and any of the parameters used as parameters of the temperature regimes (Pearson's correlation coefficients, r: -0.133 to 0.234, p>0.1). Similarly, in order to test the putative effect of the aphid biomass on the log-transformed relative proportion of the symbionts (groEL copy number / atpD copy number), we used an ANCOVA model with host species as a fixed factor with two levels (C. tujafilina population, Ct; C. cedri population, Cc) and the aphid biomass as a covariate. In the model, we assumed that the regression coefficient relating the dependent variable and the covariate varied with the factor level (i.e., interactions: host x aphid biomass).
The relationship between the log-transformed relative proportion of the symbionts, on one hand, and the host species and temperature regimes, on the other, was analyzed as follows. Since high communality is expected for the temperature parameters, we performed pairwise correlation analyses for each parameter under consideration in order to select a sound number of independent variables and avoid noise introduced by redundancy. Using p≤0.01, we identified two groups of variables, one composed by T max , T min , and T ave (r: 0.735 to 0.949), and the other by T osc and SD T (r>0.96). The maximum r was 0.37 for parameters in different groups (not significant). These results keep for the two sets of measurements (C. cedri and C. tujafilina populations). We selected T ave and T osc as representatives of the first and second groups, respectively, because they have a straightforward interpretation. T ave and T osc were then used as covariates for an ANCOVA, similar to that described above, except for the covariates included in the model. In this case, we tested for the interactions host x T ave and host x T osc .
In order to examine changes causing variations in the relative proportion of the symbionts, the response of both symbiont abundances to temperature was investigated using a linear regression analysis. Four regression analyses (two hosts x two symbionts) were computed on log abundances using T ave as the independent variable.
All analyses were performed using SPSS (IBM Corp.

Results and Discussion
Aphids from the genus Cinara infest coniferous trees. In cold climate regions, most species have holocyclic lifecycles. They overwinter as eggs that hatch in the spring into asexual females, which produce a succession of several generations of viviparous parthenogenetic females during the favorable seasons. In the fall, they produce a single annual generation of sexual individuals, which mate and lay overwintering eggs. However, in warm climate regions such as the Mediterranean area, some species are anholocyclic, and the adult parthenogenetic females survive in suboptimal temperature conditions (too cold or too warm) by migrating onto the plant's trunk or to the roots, the temperatures at which are milder (35). C. cedri lives in compact colonies on the twigs, branches, and trunks of cedars. When the temperature is too cold or too warm, adults hide under the tree bark or on the roots, and reappear at the beginning of spring. C. tujafilina mainly infests the twigs of thuja trees. In studies performed in Poland, adult aphids were found in large numbers in spring, but also at the beginning of autumn when temperatures range between 15 and 25°C, with an optimum at 25°C. The nymphs migrate to the roots when adverse temperatures are main-tained for more than 72 h, while the adults remain on the trunk until death (10,11). In the present study, even though we searched for aphids on the plant stem and leaves twice a week in 2012, and weekly in 2013, we were only able to detect C. cedri between May and June in 2012 and in June, 2013 (Fig. 1A). As for C. tujafilina, we were able to collect aphid samples for a longer period of time, between May and July and between October and December in 2012, and between April and July in 2013 (Fig. 1B). An additional sample was collected in November that year. Differences in the presence of aphids between 2012 and 2013 may have been due to the unusual weather conditions in 2013 in the Mediterranean area, according to data collected by AEMET, with a cool spring, which delayed the onset of C. cedri colonies on the tree branches, and warm autumn, which delayed the migration of C. tujafilina to the twigs after summer.
We determined endosymbiont densities in our samples using qPCR with the single copy genes groEL and atpD as indicators of the amount of each endosymbiont in the consortia in both aphid species under study. The results obtained are shown in Tables 2 and 3 (for C. cedri and C. tujafilina, respectively). We observed differences in the relative abundance of B. aphidicola and S. symbiotica depending on the host (Fig. 2 and 3). In C. cedri, B. aphidicola BCc was always present at a higher proportion, and only slight variations among samples were observed (between 80.7 and 93.6%), with an average of 89.33% in 2012 and 82.99% in 2013 ( Fig.  2A and 2B). In contrast, relative amounts were highly variable in C. tujafilina (Fig. 2C and 2D; Table 3). B. aphidicola BCt oscillated between 69.80% and 41.08% in 2012, and between 73.93% and 45.30% in 2013.
Since the individual insect size differed in the different samples, we first tested whether there was a relationship between insect weight and the relative abundance of B. aphidicola and S. symbiotica in both aphid species. We found that the two aphid species populations analyzed differed significantly in their average values for the log-ratio 'Buchnera/Serratia' (ANCOVA, p<0.001) after controlling for aphid weight, whereas the log-ratio did not correlate with aphid weight in the two species under study (ANCOVA, p=0.372 for Ct, 0.057 for Cc).
We then investigated the influence of temperature on the relative amount of both endosymbionts, using the daily average temperature and temperature oscillations as independent variables. The results obtained revealed that the two aphid species populations analyzed differed significantly in their average values for the log-ratio 'Buchnera/Serratia' (ANCOVA, p<0.001) after controlling for T ave and T osc . However, while the changes observed in the log-ratio did not correlate with T osc in the two species (ANCOVA, p>0.175), the relationship of this log-ratio with T ave was dependent on the host species (ANCOVA, p=0.024; Fig. 3). In C. cedri, slight changes in this log-ratio did not correlate with T ave (ANCOVA, p=0.190), whereas these changes positively correlated in C. tujafilina (ANCOVA, p=0.016).
In the case of C. cedri, a linear regression analysis of the response of endosymbiont abundance (measured as copy numbers of groEL from B. aphidicola and atpD from S. symbiotica) to T ave indicates that the cell densities of both endosymbionts did not correlate with daily average tempera-    adapt to high temperatures, which may, in turn, result in the disappearance of this species from the Mediterranean area as a consequence of climate change.
On the other hand, the behavior of C. tujafilina indicates that it still has some potential to cope with high temperatures. Based on knowledge on A. pisum, if S. symbiotica SCt retains some protective role against heat stress, we speculate that it will tolerate the consequences of climate change in the Mediterranean area better. Burke et al. (2) analyzed pea aphid populations, with and without protective S-symbionts, exposed to heat stress. They corroborated that when A. pisum was infected with S. symbiotica, B. aphidicola cell density was not altered in response to heat exposure; however, many S. symbiotica cells simultaneously lysed. Therefore, the density of S. symbiotica was reduced to one half after a heat shock treatment at 39°C. However, in the C. tujafilina population analyzed in the present study, a statistical analysis of the absolute copy numbers of groEL and atpD revealed that although both endosymbiotic bacteria showed a trend line with a positive slope (0.0540 for B. aphidicola BCt; 0.0296 for S. symbiotica SCt; Fig. 4B), a significant positive correlation was only observed between the average temperature recorded and the amount of B. aphidicola BCt (p=0.002). This result indicates that the increase observed in the Buchnera/Serratia ratio in C. tujafilina (Fig. 3) was due to an increase in B. aphidicola BCt abundance, and not to a decrease in S. symbiotica SCt cell density. Nevertheless, the increase found in B. aphidicola BCt cell density with higher daily average temperatures may be in accordance with the hypothesis of a protective role for S. symbiotica SCt. It is important to note that, as in previous studies (2, 38), we used single copy genes as indicators of the amount of each endosymbiont in the consortia. Although the ploidy of S. symbiotica has not yet been investigated, B. aphidicola from A. pisum has the ability to possess more than 100 genome copies cell -1 (20), and the genome copy number varies in response to the developmental stage and morph of the host (21). Even though we used aphids from the same developmental stage, the possibility of an increase in the B. aphidicola BCt gene copy number due to genome amplification in response to heat stress cannot be ruled out. In any case, an increase in the number of genomes may also be linked to an increase in gene function.
Based on metabolomic analyses, it was hypothesized that the lysis of S. symbiotica after heat exposure may deliver protective metabolites that enhance B. aphidicola survival (2). The authors suggested that molecules such as N-acetyl-D-mannosamine, mannose-6-phosphate, and β-alanine are involved in this function. Since no decrease occurred in the number of S. symbiotica cells in C. tujafilina, cell lysis does not appear to release putatively protective metabolites. Nevertheless, it is important to note that most enzymes involved in the metabolism of these compounds, which are encoded in the S. symbiotica SAp genome, are also retained in strain SCt, while they are completely absent in SCc (4,23,29). However, panD, encoding aspartate 1-decarboxylase (EC 4.1.1.11), the protein involved in the biosynthesis of β-alanine from L-aspartate, has not been annotated in S. symbiotica SCt. Therefore, it is possible that the genome reduction undergone by this strain after its recent establishment as an obligate endosymbiont is already affecting its ability to protect the symbiotic consortium from heat stress. In any case, our results are only indirect evidence, and further studies using controlled experiments are needed in order to test the hypothesis.
In summary, we have performed a study on two populations belonging to two highly related conifer aphid species that present a Buchnera/Serratia consortium, but differ in the stage of their obligate relationship with the host. Our results indicate that S. symbiotica SCc, with a well-established co-obligate endosymbiotic status, has lost its ability to protect the aphid and B. aphidicola from heat stress. On the other hand S. symbiotica SCt, with an early stage of integration in its consortium with B. aphidicola, may still retain some protective role; however, the mechanisms responsible for this protection do not appear to involve cell lysis, as it was proposed in pea aphids. The difference between both consortia may explain the two aphid population dynamics through the annual seasons and may also make C. cedri more vulnerable to the increase expected in the frequency of extreme heat events in the Mediterranean area predicted by climate change models.