There is a considerable body of research on key Mexican domesticated species like maize (e.g., Arteaga et al., 2016; Romero Navarro et al., 2017), but limited work for most of the rest (Bellon et al., 2009; Piñero et al., 2009; **Figure 1**). Similarly, maize landraces' distribution and domestication history has been widely analyzed (Kato et al., 2009; Comisión Nacional para el Conocimiento y Uso de la Biodiversidad [CONABIO], 2011), but this information remains unknown or has not been systematized for the rest of the species. This initiative aims to study, conserve, evaluate, safe keep, analyze, and sustainably use this genetic diversity.

# Functional Genomics

From studies in functional genomics we know that the phenotypes obtained during the early stages of domestication are governed by relatively few genomic loci, which, in general, are different from the loci involved in later phenotypic diversification, and from the loci subjected to natural selection (Meyer and Purugganan, 2013). Research on these topics has mostly been conducted abroad, but Mexican research institutions devoted to these areas were recently created. In developing this research, it is important to consider native varieties and molecular markers representative of this variation (Ganal et al., 2011; Caldu-Primo et al., 2017).

This component aims to have well-annotated genomes and diversity panels for the main agriculture and forestry species of importance to Mexico, along with their wild relatives. This should allow for a better understanding of the molecular basis of domestication, but also to apply this information to breeding, conserving, and monitoring.

# Production

In Mexico, most agricultural and forestry land belongs to campesino or indigenous communities. Their agriculture tends to be performed in blocks <5 ha, which characterizes them as smallholders. Although they are often perceived as 'unproductive,' the campesinos' aggregated production is the backbone of Mexico's food security (Bellon et al., under review). However, since the 1980s, the programs directed at smallholder agriculture were drastically reduced, or redirected to provide farmers with Green Revolution packages (Turrent-Fernández and Cortés-Flores, 2005b). Remarkably, many smallholders continue to use traditional varieties (Eakin et al., 2014), because these are adequate and competitive when grown under their local conditions (Muñoz et al., 1976), and because smallholders not only focus on yield, but also look to fulfill cultural preferences (Brush and Perales, 2007; Bellon and Hellin, 2011).

Importantly, smallholder farmers tend to be able to obtain some degree of usable yield in underperforming environments without aid (Bellon et al., under review). This has two main implications. First, the productivity gap of many of these farmers could be closed with minor agronomic improvements and breeding, with important consequences. Second, millions of farmers, spread across a wide variety of environments and cultural preferences, represent the ideal scenario for evolution under domestication to continue to occur at the scale and range needed to effectively maintain and generate new genetic diversity.

This component aims attending agricultural and forestry production so that they can be increased in a sustainable and

sovereign way, while also fulfilling the Mexican needs regarding desired cultivars, quantity, quality, and local cultural preferences. For this, it is crucial to protect smallholders' production and target it with adequate programs that consider management, markets, and education.

knowledge for each species. Differences in line color are only for visual purposes.

# Plant Breeding

Formal plant breeding was introduced into Mexico in replication of the model of the United States. However, in Mexico most production does not use seeds generated by 'formal' breeding (Turrent-Fernández and Cortés-Flores, 2005a,b; Donnet et al., 2012). This could be interpreted as a failure for formal breeding, but also as the success of smallholders' practices of domestication.

Commercial lines have not been widely adopted for several motives, one is that the formal breeding lines tend not to be useful to the smallholders' production systems and environments. This is a consequence of the way breeding is performed and targeted. The environmental conditions characteristic of the areas where the vast majority of smallholders are located are poorly represented by the Mexican public research stations where improved lines have been developed (Bellon et al., 2005). Also, the selection goals of the smallholders can vary depending on the particular characteristics of the local environment and the producer's interests, which do not necessarily focus on maximizing yield. For these reasons, in the areas and conditions where local varieties are traditionally grown, they tend to have a better performance than the improved lines in terms of yield, nutritional value, forage quality, local appreciated taste, or precocity (Muñoz et al., 1976; Perales et al., 1998; Sociedad Mexicana de Fitogenética [SOMEFI], 2007, 2009).

Given the range of conditions under which smallholders conduct agriculture and the limited success of formal breeding over the last 70 years in Mexico, it is unlikely that commercialtype breeding would become a successful strategy. Instead, Mexico should recognize, and incorporate into breeding, the diversity of traditional and local knowledge on shaping and adapting cultivars. For this, public policies should be shifted to better support informal breeding, and to focus formal breeding

on the smallholders' needs and adaptation to local conditions. For this, landraces should be incorporated as the base material of breeding programs, instead of only as donors for elite materials (Comisión Nacional para el Conocimiento y Uso de la Biodiversidad [CONABIO], 2017).

The aim for this component is to have breeding programs for a range of native species, using alternative tools to accelerate and improve the breeding process, where the objective would be breeding for the smallholders' needs under present and future, social and environmental conditions. To accomplish this, alternative tools need to be explored, including documenting and sharing campesino-to-campesino experiences, participatory breeding, genomic selection and evolutionary breeding.

# Environmental Change

Mexico's agriculture and forestry production are facing environmental change in the form of soil degradation, pollution, invasive species and climate change. It is difficult to generalize how environmental change would affect a particular crop or wild species; for instance, the effect of climate change on maize depends on the plant's genotype, local environment, and management (Mercer and Perales, 2010). Nonetheless, we know that environmental change has important economic impacts. For example, ∼10% of Mexican agriculture land is eroded, which translates into ∼50% of PROCAMPO aid costs (Sánchez-Colón et al., 2009; Ávalos et al., 2011). Given the large diversity of environments where Mexican crop and wild species occur, it is likely that useful variation to cope with new sources of stress already exists. What is needed is to make that diversity available to producers by enhancing seed-exchange networks, breeding and access to seeds, and environmental information at local and national scales. Therefore, this component aims to accelerate the adaptation of cultivated plants to environmental change and look for effective mechanisms to conserve the capacity of wild populations to adapt.

# ROADMAP

Research on the previous components could target and connect both the theoretical understanding of the evolutionary processes, and the practical aspects of applying this knowledge to conservation and production. However, for this to happen in a National scale, it is necessary to systematize and make data available to both the academy and wider public, and to influence public policy. CONABIO is a Mexican inter-ministerial commission in charge of that type of activities regarding biodiversity. We therefore envision that CONABIO would help to develop the initiative's in the following early stages:

Examples of the first steps are already ongoing. For instance, in the case of wild relatives, a systematic conservation planning analysis is being performed incorporating both the distribution of genetic diversity and social variables1 . For the cultivated forms, maize is the species with more data available on the distribution of native races and their genetic characterization. What is next needed is to integrate its management practices, uses and environmental and microbiome data, to then make this information available and accessible to farmers, breeders and wider audience. The target of this should be to strengthen seed exchange networks, resources for participatory breeding and campesino-to-campesino experience sharing. Molecular tools would help to accelerate breeding, especially at the stage of crosses design and genetic diversity monitoring. We estimate that with minimal breeding support and agronomic improvements, it would be possible to increase the average yield of 4 million ha from 1.3 to 2.3 ton/ha. This would be enough to cover the maize needs of c.a. 88.5 million people (Comisión Nacional para el Conocimiento y Uso de la Biodiversidad [CONABIO], 2017) without implementing intensive agriculture systems.

# AUTHOR CONTRIBUTIONS

AM-Y, MCR, and CB performed literature reviews and contributed to the discussion. AM-Y, DP, FAG, and JS wrote the manuscript. All authors conceived and designed the manuscript.

# FUNDING

This work was supported by Secretaría de Medio Ambiente y Recursos Naturales (SEMARNAT) through the grant "Contribución de la Biodiversidad para el Cambio Climático" to CONABIO and by Consejo Nacional de Ciencia y Tecnología through the grant 247730 to DP.

# ACKNOWLEDGMENTS

This document represents a summary of CONABIO's ideas and aims, which we are constructing in collaboration with external researchers, organizations, and areas within CONABIO. We are thankful to all these colleagues for helping to conceive the initiative, to Mauricio Bellon and Nancy Arizpe for

1http://www.psmesoamerica.org/en/

improving it regarding social aspects, to Patricia Koleff, Elleli Huerta, Raúl Jiménez, Rafael Obregón, and Jorge Larson for their comments to this document, to María Andrea Orjuela, Oswaldo Oliveros, and Alejandro Ponce for the technical contributions, and to Karl Philips for editing the English of the manuscript.

# REFERENCES

# SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2018.00209/ full#supplementary-material

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The handling Editor declared a shared affiliation, though no other collaboration, with the authors.

Copyright © 2018 Mastretta-Yanes, Acevedo Gasman, Burgeff, Cano Ramírez, Piñero and Sarukhán. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# The Mating System of the Wild-to-Domesticated Complex of Gossypium hirsutum L. Is Mixed

Rebeca Velázquez-López1 \* † , Ana Wegier1 \* † , Valeria Alavez1 , Javier Pérez-López1 , Valeria Vázquez-Barrios1 , Denise Arroyo-Lambaer1 , Alejandro Ponce-Mendoza2 and William E. Kunin3

1 Laboratorio de Genética de la Conservación, Jardín Botánico, Instituto de Biología, Universidad Nacional Autónoma de México, Mexico City, Mexico, 2 Comisión Nacional para el Conocimiento y Uso de la Biodiversidad, Mexico City, Mexico, 3 Department of Ecology and Evolution, Faculty of Biological Sciences, University of Leeds, Leeds, United Kingdom

#### Edited by:

Charles Roland Clement, National Institute of Amazonian Research, Brazil

#### Reviewed by:

Sevan Suni, Harvard University, United States Alexandre Magno Sebbenn, Instituto Florestal, Brazil Fernanda Amato Gaiotto, Universidade Estadual de Santa Cruz, Brazil

#### \*Correspondence:

Rebeca Velázquez-López<EMAIL_ADDRESS>Ana Wegier<EMAIL_ADDRESS><EMAIL_ADDRESS>
†These authors have contributed equally to this work.

#### Specialty section:

This article was submitted to Evolutionary and Population Genetics, a section of the journal Frontiers in Plant Science

> Received: 18 August 2017 Accepted: 13 April 2018 Published: 09 May 2018

#### Citation:

Velázquez-López R, Wegier A, Alavez V, Pérez -López J, Vázquez -Barrios V, Arroyo-Lambaer D, Ponce-Mendoza A and Kunin WE (2018) The Mating System of the Wild-to-Domesticated Complex of Gossypium hirsutum L. Is Mixed. Front. Plant Sci. 9:574. doi: 10.3389/fpls.2018.00574 The domestication syndrome of many plants includes changes in their mating systems. The evolution of the latter is shaped by ecological and genetic factors that are particular to an area. Thus, the reproductive biology of wild relatives must be studied in their natural distribution to understand the mating system of a crop species as a whole. Gossypium hirsutum (upland cotton) includes both domesticated varieties and wild populations of the same species. Most studies on mating systems describe cultivated cotton as self-pollinated, while studies on pollen dispersal report outcrossing; however, the mating system of upland cotton has not been described as mixed and little is known about its wild relatives. In this study we selected two wild metapopulations for comparison with domesticated plants and one metapopulation with evidence of recent gene flow between wild relatives and the crop to evaluate the mating system of cotton's wild-to-domesticated complex. Using classic reproductive biology methods, our data demonstrate that upland cotton presents a mixed mating system throughout the complex. Given cotton's capacity for outcrossing, differences caused by the domestication process in cultivated individuals can have consequences for its wild relatives. This characterization of the diversity of the wild relatives in their natural distribution, as well as their interactions with the crop, will be useful to design and implement adequate strategies for conservation and biosecurity.

Keywords: cotton, crop wild relatives, mating system, reproductive success, xenogamy, autogamy, domestication process, introgression

# INTRODUCTION

Plant domestication is a complex and continuing process (Casas et al., 2007; Vaughan et al., 2007). For 10,000 years, humans have selected attributes of interest in a range of economically valuable plants through their management and utilization (Gepts, 2004); consequently, different techniques, trait preferences, environments, and selection intensities have shaped the degree of domestication of each species (Meyer and Purugganan, 2013). Today we can find: (1) crop populations that are highly domesticated and depend on human intervention for survival; (2) semi-domesticated populations with recognizable traits of the domestication syndrome, but able to survive in the wild if human intervention ceases; (3) incipiently domesticated populations whose selected traits

**269**

have not yet diverged markedly from those found in wild populations; (4) incidentally co-evolved populations that adapt to human disturbed environments, but without direct human selection; (5) feral populations derived from 2, 3, or 4; and (6) wild relatives (Clement, 1999). Given these diverse scenarios, the biological diversity contained in wild-to-domesticated complexes should be considered in studies about crop ecology and evolution (Warwick and Stewart, 2005; Casas et al., 2007).

In plants, one of the key life history traits is the mating system (Vaughan et al., 2007). This feature helps determine the genetic composition of populations and, therefore, has a crucial role in the evolution of species (Charlesworth, 2006); additionally, it explains who is mating with whom, which is a fundamental issue for conservation biology (Barrett and Harder, 1996). The mating system often changes during domestication (Meyer et al., 2012), and wild relatives contain the plesiomorphic state of this trait (Ellstrand et al., 1999; Doebley et al., 2006; Andersson and de Vicente, 2010). A shift from the ancestral system toward a new one can be selected until fixation; for instance, there are some crops that are unable to reproduce without human intervention (Ellstrand et al., 1999), such as vegetatively propagated sycamore fig and other fruit trees (Zohary and Spiegel-Roy, 1975). Some crops have multiple mating systems, such as domesticated 'Maradol' (Carica papaya), which is hermaphroditic, while native varieties and wild papayas are dioecious (Carvalho and Renner, 2012). Importantly, the characterization of the mating system of many plant species has been biased toward the domesticated counterparts, because only a sub-sample of the wild-to-domesticated complex was used (e.g., Carica papaya (Damasceno et al., 2009), Persea americana (Ish-Am et al., 1999), Piper nigrum (Thangaselvabal et al., 2008). This bias may have profound consequences for the conservation of conspecific wild relatives, especially because conclusions drawn from studies with domesticated varieties are extrapolated to the whole species, failing to consider the genotypic and phenotypic diversity that wild relatives possess. The conservation of this diversity is fundamental, because it is a genetic reservoir that includes a wider range of adaptive traits that may be of additional agricultural relevance, such as resistance to pests and pathogens and tolerance to abiotic stresses (Warschefsky et al., 2014).

Upland cotton, Gossypium hirsutum, is an economically important plant species, particularly known for being the leading source of natural fiber. Worldwide, over 90% of cotton production comes from cultivars of G. hirsutum and in 2014 the species ranked eighth in the world's harvested area, reaching almost 35 million hectares (Crop production, FAOSTAT, 2017). Given the economic importance of the species, its mating system has been the focus of several studies since 1903 (Simpson, 1954); however, the majority of them concentrated on domesticated cotton and have described it as predominantly autogamous and self-pollinated (see Supplementary Material 1). On the other hand, studies on pollen dispersal of G. hirsutum, from the beginning of its modern breeding as a crop to the present day, refer to cotton's ability to produce offspring by crossing (1944– 2016; see Supplementary Material 1 for a review); however, the mating system is not described as mixed (Loden and Richmond, 1951; Richmond, 1951; Simpson, 1954; Imam and Allard, 1965; Meredith and Bridge, 1973). A specific study on the mating system of wild populations in their natural distribution is lacking.

In Mesoamerica, G. hirsutum exists as a complex of wild to domesticated forms (Brubaker and Wendel, 1994); hence, it is an ideal region to characterize the mating system of this upland cotton complex, identify possible differences, and to integrate this information into regional management plans. In Mexico - its center of origin, diversity and domestication (Ulloa et al., 2005; Burgeff et al., 2014; Pérez- Mendoza et al., 2016) – the complex includes cultivated and highly improved varieties, genetically modified varieties, traditionally managed landraces, feral, and wild populations. All of them belong to the primary gene pool of the species (Andersson and de Vicente, 2010) and gene flow among them occurs, even over long distances (Wegier et al., 2011). Moreover, eight wild G. hirsutum metapopulations have been recognized, based on geographic, ecologic, and genetic differences (Wegier et al., 2011; Bauer-Panskus et al., 2013). Wegier et al. (2011) demonstrated that recent gene flow, followed by introgressive hybridization, occurs between a number of wild populations distributed in the north and south of Mexico, and commercial cotton cultivars in the northern states of the country. Our study provides the first data on the mating system of wild G. hirsutum in situ within the natural distribution of the species in Mexico. In order to assess this, we evaluated the capacity of domesticated cotton, wild cotton, and wild cotton with evidence of introgression, to produce offspring by either xenogamy (cross-pollination between different genets) or autogamy (self-fertilization).

# MATERIALS AND METHODS

# Study System

Upland cotton, G. hirsutum L., is a species with wild, feral, and semi-domesticated populations (Brubaker and Wendel, 1994). All cultivated forms, including the highly improved varieties or genetically engineered varieties, cannot be considered as fully domesticated, because they are able to survive even if human intervention stops.

Wild G. hirsutum flowers all year round. Flowers are white, hermaphrodite, cup shaped, with a single central style surrounded at the bottom by stamens (Meade, 1918; Smith and Cothren, 1999). Some plants exhibit flowers with a colored disk inside of the base of the cup that ranges from deep red to light yellow (Tan et al., 2013). Flowers remain open between 8 and 11 h; at the start of the day, they are all white and when they close the sepals start turning pink at the base (Smith and Cothren, 1999). Anthesis takes place in the morning, as soon as the flower completely opens, and the stamens start to release pollen soon afterward (Smith and Cothren, 1999). Flowers produce both pollen and nectar as a reward for visitors (Wäckers and Bonifay, 2004).

# Study Sites

Sampling was performed in coastal dunes and dry forests of Mexico, in three of the eight wild cotton metapopulations defined genetically, geographically, and ecologically by

Wegier et al. (2011), namely: Central Pacific Metapopulation (CPM), Yucatan Peninsula Metapopulation (YPM), and South Pacific Metapopulation (SPM). SPM is of particular interest due to evidence of recent introgressive hybridization with domesticated plants (Wegier et al., 2011). Given the distinctive extinction-colonization dynamic observed in cotton metapopulations (Wegier, 2013), the full extent of CPM, YPM, and SPM was surveyed to find G. hirsutum patches with enough flowers (**Figure 1**). Hand-pollination treatments (Tate and Simpson, 2004; Machado and Sazima, 2008; Hernández-Montero and Sosa, 2016) were carried out during the dry season, between November 2012 and May 2013: CPM in November 2012, YPM in December 2012, and SPM in February 2013. Sites were revisited for fruit collection after 86 days, on average. In addition, domesticated cotton plants, bought in local markets, were kept under greenhouse conditions in Mexico City to maintain a suitable temperature (**Figure 1**).

# Mating System

In order to execute the hand-pollination treatments to test for different mating systems, a search was conducted for flower buds before anthesis (Tate and Simpson, 2004; Machado and Sazima, 2008; Hernández-Montero and Sosa, 2016). In each metapopulation, 40 replicates of the five pollination treatments were set up (**Table 1**), anticipating the risk of collecting too few fruits afterward: assisted self-pollination, automatic selfpollination, assisted cross-pollination (Kearns and Inouye, 1993), emasculated control (to avoid automatic self-pollination), and control (open-pollination). Multiple treatments were placed on the same plant where possible to control for individual variation; however, due to the variability of the number of flowers, not all of the plants held the same type or number of treatments. Moreover, when flowers were scarce, treatments were placed daily within each study site, in up to four patches per metapopulation, until the 40 replicates per treatment were completed. Special care was taken to avoid changing or altering the environment (i.e., without introducing new genotypes or changing plant abundances or distributions). The same experimental design was applied for domesticated plants in a greenhouse that allowed the entry of local insects. Some treatments required mesh bags to exclude any pollinator access that could alter the results (**Table 1**). The treatments that did not include bagging before anthesis were bagged after flower closure to help control for mechanical damage from bagging.

# Reproductive Success

Fruit-set was calculated as the percentage of recorded fruits produced by each treatment in each metapopulation (Dafni, 1992). In each study site, 20 flowers not involved in the pollination-treatments were collected and brought back to the laboratory in separate sealed containers with 70% alcohol. Each flower was dissected to count the number of ovules present. An average number of ovules was calculated for each wild metapopulation and for domesticated plants. Afterward, seed-set (Schoper et al., 1987; Burd, 1994) was calculated as the percentage of seeds obtained from each fruit for each pollination treatment in relation to the average number of ovules of the study population to which these fruits belonged. Additionally, all seeds were weighed individually to estimate the seed weight per treatment in each study site. Later, all the seeds were germinated individually. Each seed was washed with 2% Captan (PESTANAL R , Merck) solution and covered with a damp cotton swab; tissue culture lids were used. Seeds were checked daily until all reached emergence of the radicle. While some studies on seed germination consider only a set of seeds (Schemske, 1983; Gil and López, 2015; Raphael et al., 2017; Farooq et al., 2018), we took into account all of the collected seeds for the analysis.

# Outcrossing Rate

The outcrossing rate (Te) was calculated for each study site following Barrett et al. (1996):

$$Te = 1 - S \tag{1}$$

where S is the selfing rate, estimated with the fruit-set results from our selfing (Ws) and outcrossing (Wx) treatments, i.e., automatic self-pollination and emasculated control, respectively. For CPM, Wx was obtained with the fruit-set from the assisted crosspollination treatment, because none of the emasculated control results were found when revisiting the metapopulation for fruit collection.

$$\mathbf{S} = (\frac{\boldsymbol{\omega}\_{\mathbf{s}}}{\boldsymbol{\omega}\_{\mathbf{s}} + \boldsymbol{\omega}\_{\mathbf{x}}}) \tag{2}$$

# Statistical Analyses

To test if there were significant differences in seed-set and seed weight among treatments, a Generalized Linear Mixed Model GLMM (Zuur et al., 2009) was used considering the plant as a random factor, because the pollination treatments were not equally represented in each plant (as explained in section 2.3). For GLMM analyses, a Quasi-Poisson distribution was considered for seed-set and a Gaussian distribution for seed weight (Cayuela, 2009). Afterward, a Tukey post hoc test was performed to evaluate the significance of the results. To compare germination frequencies and percentage of fruit-set, a chi square test was used with the post hoc standardized residue test for each one. Outliers were identified using the method described by Viechtbauer and Cheung (2010); to summarize, a multivariate detection method (Cook distance) was used to calculate the distance among all data points, and those that were not included in the general model were identified as "influential data points" or outlier values. Germination was calculated as the number of germinated seeds in relation to the total number of seeds (Gómez, 2004; Gil and López, 2015). All tests were carried out with the lme4, multcomp, stats, and ggplot2 packages of R version 3.4.3 (R Core Team, 2017). The scripts utilized for the analyses are available online at https://github.com/conservationgenetics/ BiologiaReproductiva.git.

# RESULTS

# Fruit-Set, Seed-Set, and Seed Weight

All treatments produced fruits regardless of the metapopulation (**Table 2**). The CPM open-pollinated control showed the highest

FIGURE 1 | Approximate distributions of the Central Pacific (blue), Southern Pacific (green) and Yucatan Peninsula (orange) upland cotton metapopulations studied in Mexico, with locations of field study sites (leaves) and the greenhouse in Mexico City where cultivated upland cotton was grown out for this study. The Central Pacific and Yucatan Peninsula metapopulations are wild and the Southern Pacific metapopulation contains populations with evidence of recent introgression with cultivated upland cotton.

TABLE 1 | Characteristics of each of the pollination treatments applied to G. hirsutum flowers, following Dafni (1992), in three field wild metapopulations and cultivated cotton in the greenhouse.

value of all treatments among all groups; YPM showed the highest fruit-set produced by outcrossing, and the lowest by automatic self-pollination and control treatment. On the other hand, the highest fruit-set was observed for all the treatments in SPM, with exception of the open-pollination control.

Seeds were produced both by outcrossing and selfing treatments in all metapopulations (**Figure 2**). The average number of seeds per fruit was 15.9 in SPM, 15.4 in domesticated, 12.3 in CPM, and 10.0 in YPM, while the average number of ovules was 16.3 in SPM, 28.6 in domesticated, 15.6 in CPM, and 13.4 in YPM. Regarding seed-set, the control treatment of domesticated cotton was lower than that of wild and introgressed plants [P(χ 2 ) = 4.13 × 10−4 , df = 3] (**Figure 3**). When evaluating the results of each metapopulation individually, CPM presented seed-set differences between the control and the rest of the treatments (P(χ 2 ) = 0.1 × 10−4 , df = 3), in SPM the differences were found between cross-pollination and all treatments, except emasculated control [P(χ 2 ) = 0.001, df = 4], while in YPM and the domesticated there were no significant differences among treatments (**Figure 2**). On the other hand, seed weight only presented differences in SPM, between cross-pollination and both assisted and automatic self-pollination [P(χ 2 ) = 0.026,

∗Categories significantly different using adjusted standardized residuals greater than 2.0 and less than −2.0 (χ 2 test). ND: Fruit-set could not be determined because none of the treated flowers were found when revisiting the metapopulation for fruit collection. NA: the chi-square test could not be determined because one of the treatments was not collected.

df = 4] (**Figure 2**). In YPM, CPM, and the domesticated cotton, no significant differences among treatments were observed (**Figure 2**).

# Germination

Less than 20% of the seeds from wild metapopulations CPM and YPM germinated. Regarding the domesticated group, 40– 63% of the seeds germinated, except for the seeds produced by cross-pollination that only reached 28%. The seeds from the five treatments assessed at SPM showed germination percentages above 83%, with cross-pollination reaching the highest value of 96% (**Table 3**).

With regard to the germination rate of the seeds that germinated (**Figure 4**), the slope of the curve suggests that wild upland cotton presents some kind of inhibition to the completion of germination, whereas domesticated populations do not display this behavior. As shown in **Figure 4A**, domesticated seeds germinated faster, within the first 6 days, whereas seeds presenting evidence of introgressive hybridization (SPM) reached 95% of germination within the first 7 days and continued germinating for 48 days. Unlike domesticated and SPM seeds, the seeds of wild plants germinated over the course of 73 days (**Figure 4A**). Concerning the pollination treatments from all study sites, 50% of the seeds of all treatments germinated within the first 5 days; however, after the 5th day the difference in germination rate is evident between autogamy and the rest of the treatments (**Figure 4B**).

## Outcrossing Rate

All study groups presented outcrossing rates different from 0 and 1 (i.e., 1 > Te > 0), which is indicative of a mixed mating system. Wild metapopulations (YPM and CPM) recorded a higher outcrossing rate (0.72 and 0.71, respectively) than SPM (0.40) and domesticated (0.65) (Supplementary Material 1).

# DISCUSSION

# Mating System of Upland Cotton's Wild-Domesticated Complex

Richards (1997) defined autogamy as within-flower or selfpollination, and allogamy as the pollination between pollen and ovules of different flowers; moreover, he further divided allogamy into geitonogamy (i.e., pollination between different flowers on the same genet) and xenogamy (i.e., pollination between pollen and ovules of different genets). Our results show that wild and domesticated cotton produce offspring in all pollination treatments (**Figures 2**, **4B** and **Tables 2**, **3**); thus, the analyzed plants have the capacity to produce progeny by both autogamy and xenogamy. To discard geitonogamy, it is necessary to perform a molecular genetic analysis of paternity. However, since autogamy is common in our system, there is no need to discard this type of allogamy. Furthermore, previous studies (see Supplementary Material 1), together with our own, indicate that the G. hirsutum wild-domesticated complex has a mixed mating system. This result is particularly relevant in upland cotton's center of origin, because of its significance on strategies for long-term conservation of genetic diversity in the event of gene flow between wild and domesticated relatives (Ellstrand, 1992).

Barrett and Eckert (1990) and Barrett et al. (1996) described the outcrossing rate (Te), which indicates that when the value is 0.5 the mating system is equally balanced between self and cross-pollination. Any value different from 0 (completely selfpollinated) or 1 (completely cross-pollinated) implies a mixed mating system; when Te > 0.5, the system is predominantly allogamous-xenogamous, whereas when Te < 0.5, the system is predominantly autogamous. Our observed rates vary from Te > 0.5, e.g., 0.71 (CPM), 0.72 (YPM) and 0.65 (domesticated), to Te < 0.5, e.g., 0.40 (SPM). Domesticated plants, and wild CPM and YPM, have a greater contribution of seeds from crosspollination in the next generation, although the contribution of self-pollination is high and important, and it contributed to maintenance of genetic structure. The high contribution of selfpollinated seeds in SPM is striking, far from being similar or intermediate between wild and domesticated; local factors may be affecting the result and should be addressed in a future study.

To further explore the mixed mating system of the species, we compared the germination rate of seeds produced by different pollination treatments. We found that within the first 5 days the seeds for all treatments reach 50% of germination. After the 10th day, a notorious difference on germination rate (<15%) among autogamy and the other treatments is observed. Such discrepancy is due to the difference in number of seeds produced in each treatment (**Figure 4B**). As suggested theoretically, when germination does not differ among treatments, self-pollination

is not the cause of inbreeding depression (Charlesworth and Willis, 2009). The mating system described in our study coincides with Baker's law of reproductive assurance (Pannell and Spencer, 1998), where species that migrate long distances colonize or recolonize patches initially by self-fertilization; then, because of its perennial nature, generations overlap in the same area and plants are pollinated by close relatives or by themselves in the absence of pollinators (Kalisz et al., 2004). The information described here, agrees with the ecological and genetic evidence that describes the metapopulation dynamics of G. hirsutum, along with the ability to migrate long distances, historically and currently (Wegier et al., 2011).

In addition, Wegier et al. (2011) reported high values of gene flow among metapopulations in the same study area, which could homogenize genetic variation, but their data exhibit population structure (k = 8) and high FST. Self-pollination and cross-pollination seem to maintain the genetic diversity of the species in the wild, although crossings with domesticated members of the complex (Wendel et al., 1992; Wegier et al., 2011) or even domesticated plants of Gossypium barbadense (Brubaker et al., 1993; Brubaker and Wendel, 1994; Ellstrand et al., 1999; Ellstrand, 2014) might be contributing to these results. In addition, gene flow with feral cotton can also take place (Rache Cardenal et al., 2013; de Menezes et al., 2015).

Finally, one of the fitness components measured in plants is seed weight (Primack and Kang, 1989), due to the fact that larger seeds perform better because of the higher amount of resources they possess (Armstrong and Westoby, 1993; Westoby et al., 1996). In our research, seed weight showed no significant difference between treatments within metapopulations (**Figure 2**).

# Differences of Reproductive Traits Within G. hirsutum's Wild to Domesticated Complex

Our analyses show differences in characters linked to some of the reproductive structures of upland cotton, which can be associated with the domestication syndrome and will be discussed below.

### Ovule Number

There are significant differences in ovule number [P(χ 2 ) = 0.001, df = 2], which was initially estimated to obtain the seed-set in each population (Supplementary Material 2). Wild plants produce on average 14.5 ovules per flower, while cultivated plants produce twice as many. Several authors have described a change in ovule number as a consequence of evolutionary processes. For instance, Pasquet (1998) found that ovule number supports the physiological division of cultivated cowpeas [Vigna unguiculata (L.) Walp.] into two different groups: cultivars able to flower early under inductive conditions, with ovule number lower than 17 (Biflora and Melanophthalmus) and cultivars not able to do so, with ovule number higher than 17 (Unguiculata and Sesquipedalis). Moreover, Andargie et al. (2014) reported a pair of quantitative trait loci (QTLs; qon1 and qon3) that regulate ovule number in cowpea; the alleles from the wild parent increase this trait as opposed to the cultivated, which reveals a feature of cowpea's domestication syndrome. In the case of climbing common bean (Phaseolus vulgaris L.), among the changes that occurred during the domestication process is the modification on the number of ovules, which changed from 5–8 to 2–9 ovules (Gepts and Debouck, 1991).

## Seed-Set

As shown in **Figure 3**, there are significant differences in the seed-set of wild and domesticated populations. From a much larger number of ovules, domesticated plants (open pollination controls) produce, proportionally, a lower quantity of seeds, which implies that they are not efficiently using the resources invested on ovule production (Cilas et al., 2010). Variation in seed number per boll is produced by the interplay of the plant genetics and the environment, which in turn generates either the lack of seed fertilization or completion of embryo

TABLE 3 | Germination of seeds obtained from each upland cotton metapopulation and treatment.

+All the repetitions of the emasculated control were placed, but none of them were found when revisiting the metapopulation for fruit collection. ∗Categories significantly different using adjusted standardized residuals greater than 2.0 and less than −2.0 (χ 2 test).

growth post-fertilization (Davidonis et al., 1996); therefore, our results are influenced by the experimental design and, in the future, a common garden experiment will provide insight into the effect of the environment. In comparison, wild plants are more efficient, producing seeds from nearly all of their ovules, although the net number of seeds is smaller than that produced by domesticated fruits. Many features associated with domestication are not advantageous in terms of reproduction and survival of following generations lacking human intervention (Gepts, 2004), because the selective pressures by which they have evolved are determined by humans (see categories 1–4 of the classification proposed by Clement, 1999). As a result, gene flow between wild relatives and cultivated plants could have negative consequences (Andersson and de Vicente, 2010), however, it could also give rise to in situ reservoirs of domesticated genes for the future (Ellstrand, 2018). Each domesticated cotton plant develops 50% more descendant plants than the wild plants do within their natural distribution, so the ecological-evolutionary consequences of this result will depend on the evolutionary process and the agro-ecological or ecosystem context in which plants are developed.

### Germination

One of the traits selected for during domestication is rapid germination (Frary and Doganlar, 2003), as this helps crops to start to grow at the same time and contributes to synchronous fruiting. Over time, this trait contributes to harvesting efforts and, therefore, unconsciously selects for loss of dormancy. In natural habitats, conditions are less predictable, and dormancy will contribute to different seeds germinating in different environmental conditions (Long et al., 2015). Our results on seeds that reached germination agree with what has been described for other domesticated plants that have undergone similar evolutionary processes (Fuller and Allaby, 2009; Abbo et al., 2014; Hernández et al., 2017): domesticated seeds germinate faster and practically simultaneously, whereas their wild relatives display dormancy (**Figure 4A**).

### Distinctive Traits of SPM

With respect to SPM (selected for study because of evidence of recent introgression with domesticated plants; Wegier et al., 2011), the reproductive system is mixed, as it is in wild populations without introgression and in domesticated populations. However, some of the traits that determine reproductive success are unique to this population: the variability in seed-set values is markedly different (**Figures 2**, **3**); its fruits produce more seeds than the other populations (similar to domesticated fruits, but from half the number of ovules, which makes them very efficient) (Supplementary Material 2); and these seeds have a higher percentage of germination than the other populations (**Table 3**). These characteristics can have demographic consequences in the short term, unless there are other factors that regulate this growth. On the other hand, contrary to what was expected for SPM, their resemblance to domesticated seed germination is higher than with the wild ones. The behavior is also dissimilar for introgressed seeds, which took longer to complete germination than domesticated and wild seeds (**Figure 4A**). This last phase displays a very slow response in SPM, probably associated with the loss of physiological responses, resembling domesticated plants. It is important that these analyses are repeated in subsequent years, to confirm if there is an eco-genetic trend (Price and Waser, 1979;

Edmands, 2007; Ellstrand and Rieseberg, 2016), or if it was the result of local conditions.

# Conservation and Biosecurity Implications

Many nations want to defend the rights of the next generations to enjoy and decide about biodiversity and its services, aware that decisions made today will have an impact on the natural resources available in the future (Bennett et al., 2015; Steffen et al., 2015; Morales et al., 2017). Upland cotton is a remarkably important plant for humanity, not only due to the versatile uses of its fiber, but for many other applications (Wegier et al., 2016). It follows that cotton's wild-to-domesticated complex and its environment should be a conservation priority. Mesoamerican dry forests and coastal dunes contain the ecosystems and evolutionary processes that originated, mold, and maintain wild cotton diversity and its interactions. These evolutionary services (Faith et al., 2010; Bailey, 2011; Rudman et al., 2017) are essential for species conservation, because preserving this genetic diversity allows the capacity to adapt to environmental changes (Ellstrand, 1992; Hartl, 2000). However, the factors that mold each part of the wild-to-domesticated complex are different; for example, the conservation of native traditional varieties depends to a great extent on the communities that cultivate them, their management techniques, and interests (Zhang et al., 2007). Hence, the parts of the complex that could be used for crop improvement will depend on the objectives of the new processes of domestication and breeding (Ellstrand, 2018; Mastretta-Yanes et al., 2018).

Gene flow between crops and wild relatives should be examined on a case-by-case basis (Stewart et al., 2003), especially when genetically modified organisms (GMO) are involved, because the consequences depend on the nature of the transferred genes and their regulatory mechanisms (Ellstrand, 2003). For instance, a recent study has demonstrated that genetic modifications can affect fitness traits in the long-term (Hernández-Terán et al., 2017). An important issue to keep in mind is that for gene flow to occur, a crop must be within pollination distance of a compatible population (Ellstrand and Hoffman, 1990), but in the case of domesticated plants the distances can be shortened by human activities (Dyer et al., 2009; Wegier et al., 2011). Several studies have documented hybridization events between crops and their wild relatives; for instance, in the United Kingdom, one-third of the 36 species analyzed by Raybould and Gray (1993) hybridize with at least one element of the local flora; in the Netherlands, a quarter of 42 species does (de Vries et al., 1992); and all but one of the 13 crops reviewed by Ellstrand et al. (1999) hybridize naturally with their wild relatives in some part of their agricultural distribution (including G. hirsutum and other species of subgenus Karpas). These hybridization events could lead to a decline in wild genetic diversity, as opposed to native semi-domesticated varieties in traditional Mesoamerican systems where there is evidence that domesticated genomes have formed not only by selection under domestication, but also by gene flow with other closely related populations and species (Rendón-Anaya et al., 2017). For this reason, the wild-to-domesticated dynamics in terms of genetic diversity, reproductive biology, and gene flow should be well understood in the natural distribution of the species of interest, because extrapolating conclusions based on external or incomplete information about species complexes is inconsistent with the objectives of conservation and biosafety (Beebe et al., 1997; Acevedo et al., 2016).

In this study we found that the reproductive capacity of introgressed cotton is greater than that of wild and domesticated plants. This reveals a scenario that de Wet (1968), de Wet and Harlan (1975) and Keeler et al. (1996) had already described, where wild relatives of some introgressed crops can become weeds that are difficult to control. The wide genetic diversity of G. hirsutum, along with factors modified by traditional genetic improvement and modern genetic engineering, will be problematic for agroecosystems (Altieri, 2000) and ecosystem conservation if they increase cotton's weediness or invasiveness (Schafer et al., 2011). Cotton has already been reported to persist in a few tropical regions, such as the north of Australia, Vietnam, México, the continental United States and Hawaii (Hawkins et al., 2005; Andersson and de Vicente, 2010; USDA, 2018), so it will be necessary to monitor these changes in wild populations given the species great capacity for long distance migration by natural and anthropogenic means.

Finally, local conditions can influence the results of reproductive biology studies (Ellstrand and Foster, 1983; Hucl, 1996; Murray et al., 2002); hence, it was essential to assess the mating system of G. hirsutum within its natural distribution. Some of the factors that have an effect on the results can be associated with the environment (pollen viability, nectar production, and pollinator activity due to environmental conditions; Ahrent and Caviness, 1994; Ibarra-Pérez et al., 1997; Chaves-Barrantes et al., 2014), ecological interactions (foraging rate, floral consistency, efficiency of pollen deposition, interactions with arthropodofauna, and composition of pollinator species; Rudgers, 2004; Kessler et al., 2012; Johnson et al., 2015), as well as the landscape (species abundance and surrounding species distributions; Murray et al., 2002). In this study, the results of automatic self-pollination and the emasculated control provide evidence that autogamy and allogamy occur naturally in upland cotton's natural distribution. The occurrence of the latter highlights the importance of native pollinators on the reproductive biology of G. hirsutum and, consequently, conservation strategies should take this key interaction into consideration.

# CONCLUSION

This study found that upland cotton's wild-to-domesticated complex presents a mixed mating system. This information is new for wild, domesticated, and introgressed G. hirsutum in its natural distribution, but it is in agreement with previous studies in populations of domesticated cotton (Supplementary Material 3). Consequently, G. hirsutum should be considered as having a mixed reproductive strategy throughout its whole complex, rather than being primarily autogamous. Management strategies and policies meant to conserve the diversity of cotton's wild-to-domesticated complex must take this into account.

Furthermore, physiological differences were found between cultivated cotton and its wild relatives, especially in traits such as the number of ovules per flower, number of viable seeds per fruit, and their germination behavior. Given the evidence of gene flow and introgression, these traits should be monitored systematically in wild populations and agroecosystems of interest for conservation, as well as the impact on ecological interactions, such as pollination. On the other hand, the diversity contained in the wild-to-domesticated complex must be included in longterm conservation strategies, so that future generations can have access to genetic resources with greater chances of surviving the changing environments.

# AUTHOR CONTRIBUTIONS

AW, VA, and RV-L designed the research, participated in fieldwork, performed the analyses, and wrote the manuscript. AW coordinated the study. WK designed the fieldwork and revised the analyses. AP-M participated in fieldwork and performed the analyses. JP-L, VV-B, and DA-L conducted the analyses. All authors analyzed the results and wrote the manuscript.

# FUNDING

This work was financially supported by the project "Program for the conservation of wild populations of Gossypium hirsutum

in Mexico", DGAP003/WN003/18 funded by the DGSPRNR (Dirección General del Sector Primario y Recursos Naturales Renovables) that belongs to SEMARNAT and CONABIO; complemented with the support of CONACYT scholarships (213557, 609346) and projects CONACYT-PN247672, UNAM PAPIIT No. IV200117.

# ACKNOWLEDGMENTS

We are grateful for the support of all the people in the communities where our work was done. We express our gratitude to the reviewers, whose comments and suggestions greatly improved our manuscript, and especially thank Charles R. Clement for his recommendations and advice. Also, we wish to thank Mariana Benítez, Alejandro Correa-Metrio, Ana Elena Escalante, Ulises Rosas, and Kristy Walker for their valuable comments on our work. In addition, we would like to acknowledge the support of the colleagues who assisted with sampling and laboratory work: Néstor Chavarría, José Luis Caldú,

# REFERENCES

Erick Tovar, Atsiry López, Adriana Uscanga, Cirene Gutiérrez, Luis Barba, Adriana Calderón, Lislie Solis, Florencia García-Campusano, Haven López-Sánchez, Brian Urbano, and especially Diana Peña. Pamela Rodríguez collaborated with the historical analysis. We would also like to thank Professor Zenón Cano for his teachings that have contributed to our professional development through generations. DA-L is carrying out a postdoctorate funded by CONACYT-PN247672. JP-L and VV-B thank CONACYT for grants 609346 and 477713, respectively as well as the Maestría en Ciencias Biológicas, UNAM. Finally, we are especially grateful to Nancy Corona and the CARB-CONABIO team, who have provided support and wisdom through our entire journey.

# SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2018.00574/ full#supplementary-material

of ovules per ovary, number of seeds per pod, and seed weight. Tree Genet. Genomes 6, 219–226. doi: 10.1007/s11295-009-0242-9

Zuur, A. F., Ieno, E. N., Walker, N. J., Saveliev, A. A., and Smith, G. M. (eds). (2009). "GLMM and GAMM," in Mixed Effects Models and Extensions in Ecology with R (New York, NY: Springer), 323–341.

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Velázquez-López, Wegier, Alavez, Pérez-López, Vázquez-Barrios, Arroyo-Lambaer, Ponce-Mendoza and Kunin. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.


**3. Results and discussion**

#### **3.1. Geometric changes**

The structures of {η2 -(X@Cn)}ML2 complexes for M = Pt, Pd; X = F− , 0, Li+ and n = 60, 70, 76, 84, 90 and 96 were fully optimized at the M06/LANL2DZ level of theory. The geometries that are obtained are illustrated in **Figure 1**. The key structural parameters of the stationary points are listed in **Table 1** (the structural parameters for n = 70, 76, 84, 90 and 96 are presented elsewhere). For the Pt-C60 complex in the absence of encapsulated ions, the respective lengths of the metal-carbon bonds are 2.12 Å and 2.12 Å (**Table 1**). When the Li+ ion is encapsulated into the cage, the metal-carbon bonds remain unaltered and the respective distances between

How Important is Metal-Carbon Back-Bonding for the Stability of Fullerene-Transition Metal Complexes?... http://dx.doi.org/10.5772/intechopen.70068 115

**Figure 1.** Optimized geometries for {η2 -(X@C60)}ML2 (M = Pt, Pd; X = Li+ , 0, F− ).

is the sum of the deformation energy of A (∆E(DEF)A, which is defined as the energy of A in

tion energy term, ∆E(INT)A(BC), is the interaction energy between A and (BC) for their respec-

Advanced EDA unites the natural orbitals for chemical valence (NOCV), so it is possible to separate the total orbital interactions into pairwise contributions [23]. The advanced EDA (i.e., EDA-NOCV) further divides the interaction energy (∆E(INT)) into three main components: ∆E(INT) = ∆Eelstat + ∆EPauli + ∆Eorb. It is used for a quantitative study of π back-bonding to fullerene ligands that uses the M06/TZ2P level of theory with the ADF 2016 program package [24]. The relativistic effect is considered by applying a scalar zero-order regular approximation (ZORA) [25]. The interaction energy and its decomposition terms are obtained from a singlepoint calculation using the M06/TZ2P basis set from the Gaussian 09 optimized geometry.

complexes for M = Pt, Pd; X = F−

90 and 96 were fully optimized at the M06/LANL2DZ level of theory. The geometries that are obtained are illustrated in **Figure 1**. The key structural parameters of the stationary points are listed in **Table 1** (the structural parameters for n = 70, 76, 84, 90 and 96 are presented elsewhere). For the Pt-C60 complex in the absence of encapsulated ions, the respective lengths

into the cage, the metal-carbon bonds remain unaltered and the respective distances between

of the metal-carbon bonds are 2.12 Å and 2.12 Å (**Table 1**). When the Li+

) and B (∆E(DEF)B)). The interac-

, 0, Li+

and n = 60, 70, 76, 84,

ion is encapsulated

the product relative to the optimized isolated structure (A0

@Cn)}PtL2 .

114 Fullerenes and Relative Materials - Properties and Applications

tive optimized product structures.

**Scheme 2.** Basic EDA for {η2

**3. Results and discussion**

**3.1. Geometric changes**

The structures of {η2

C1 , C2 and Li+ are 2.29 Å and 2.29 Å. As the encapsulated ion is changed to F− , the metal-carbon bonds remain almost unchanged (2.13 and 2.13 Å), but the distance between C1 , C2 and encapsulated ions (F− ) increases (3.18 and 3.18 Å). The Li+ ion is located at a site that is close to the transition metals because of electrostatic interaction. The metal-coordinated carbon atoms of C60 are negatively charged because there is π back-donation from the metal center. For the Pt-C60 complex without encapsulated ions, the natural population analysis (NPA) shows that the atomic charges on the C1 (C2 ) atoms are −0.27 (−0.27). When the cage is encapsulated by a Li+ ion, the atomic charges on the C1 (C2 ) atoms are increased to −0.32 (−0.32) and the atomic charge on the Li atom is +0.86. Therefore, the encaged Li+ ion is attracted toward these negatively charged C atoms. However, as the encapsulated ion is changed to F− , NPA shows that the atomic charges on C1 (C2 ) atoms are decreased to −0.23 (−0.23) and the atomic charge on the F atom is negative (−0.93), so the encaged F− ion is repelled by the negatively charged C atoms. In terms of Pd-C60 complexes, it is worthy of note that the geometrical distances are generally similar to the corresponding distances for Pt-C60 complexes, but the charge distributions are different. Specifically, the encaged Li atom has a charge (+0.86) but the charges on C1 (C2 ) atoms are reduced to −0.27 (−0.27). The negative charges on metal-coordinated carbon atoms are also less for X = 0 and F− . Similar geometric changes and charge distributions are seen for n = 70, 76, 84, 90 and 96 and are presented elsewhere.

is more stable. The relative thermodynamic stability increases in the order: ∆E(X = F−

(M = Pt, Pd) at M06/LANL2DZa,b

F− 51.2 (41.0, 10.2) −14.8 −67.0 +33.8 −14.7 +20.1 0 66.0 (52.1, 13.9) – −100.8 – −34.8 – Li+ 74.0 (55.0, 19.0) 8.0 −136.6 −35.9 −64.2 −29.4

F− 30.7 (25.9, 4.8) −12.6 −45.8 29.0 −13.9 +17.6 0 43.3 (35.0, 8.3) – −74.8 – −31.5 – Li+ 51.1 (37.9, 13.2) 7.8 −109.3 −34.5 −59.2 −27.7

is large, so the latter must be responsible for an increase in thermodynamic stability. Similar results are obtained for Pd-C60 complexes, but, the distortion in fragments A or B is smaller than that for the Pt-C60 complex, as is the interaction energy, so the complex is less stable.

dThe reaction energy without zero-point energy (ZPE) correction for the product, relative to the corresponding reactants.

) = −34.5 > ∆E(M = Pt, X = Li+

action between the metal fragment and the X@Cn fragment and |∆∆E(DEF)| is smaller than |∆∆E(INT)A(BC)|. Therefore, an increase in the cage size has no effect on the basic EDA results.

In an earlier study by the authors, structural parameters and spectral characteristics were

quency (Δν) and the chemical shift (Δδ) were used to describe the character of the π-complex s. In this study, the strength of the π back-bonding strength is estimated from an energetic viewpoint using an advanced EDA method. This analysis shows the effect of encapsulated

In an earlier discussion (Section 3.2), it was proven that thermodynamic stabilities increase

), as shown in **Table 2**. In addition, |∆∆E(DEF)| is small and |∆∆E(INT)A(BC)|

.

**∆∆E(DEF)c ∆E(INT)A(BC) ∆∆E(INT)A(BC)**

How Important is Metal-Carbon Back-Bonding for the Stability of Fullerene-Transition Metal Complexes?...

(X = 0) < ∆E(X = Li+

**Table 2.** Basic EDA for {η2

**X ∆E(DEF)**

Energies are given in kcal/mol.

{η2

{η2

a

b

c

**(∆E(DEF)A, ∆E(DEF)B)**

A and B respectively represent the metal fragment and C60 cage.

The difference is relative to corresponding quantity at X = 0.

= PPh3

in the order: ∆E(X = F−

For instance, ∆E(M = Pd, X = Li+

size increases (n = 70, 76, 84, 90 and 96), the encapsulated Li+

**3.3. Advanced energy decomposition analysis (advanced EDA)**

) complexes [10]. The changes in bond length (Δr/r0

) < ∆E(X = 0) < ∆E(X = Li+

used to estimate the strength of π back-bonding for {η2

ions, metal fragments and cage sizes on π back-bonding.

*3.3.1. The effect of encapsulated ions on π back-bonding*

) < ∆E

, L

) = −35.9 kcal/mol. When the cage

**c ∆Ed ∆∆Ec,d**

117

http://dx.doi.org/10.5772/intechopen.70068

ion still induces a stronger inter-

(M = Pt, Pd; X = 0, Li+

), bond angle (Δθav), vibrational fre-

), because the interaction energy (∆E(INT))

**Table 1.** Selected geometrical parameters (bond distances in Å) and the NPA atomic charge for optimized complexes ({η2 -(X@C60)}ML2 ) at the M06/LANL2DZ level of theory.

#### **3.2. Basic energy decomposition analysis (basic EDA)**

In order to better understand the factors that govern the thermodynamic stability of {η2 -(X@Cn)} ML2 complexes, basic EDA is performed on {η2 -(X@Cn)}ML2 complexes (the basic EDA results for n = 70, 76, 84, 90 and 96 are presented elsewhere). The bonding energy (∆E) is defined as ∆E= E({η2 -(X@Cn)}ML2 ) − E(X@Cn) − E(ML2 ), using Eq. (1). The Pt-C60 complex without encapsulated ions is initially considered. Basic EDA shows that both the metal fragment and the empty C60 fragment are distorted during the formation of the metal-carbon bond (**Table 2**). The metal fragment undergoes greater distortion (∆E(DEF)A = 52.1 kcal/mol) than C60 (∆E(DEF)B = 13.9 kcal/mol). The same results are obtained for X = Li+ or F− . It is found that the encapsulation of the Li+ ion induces more distortion in fragments A and B of the Pt-Li+ @C60 complex than those of the Pt-F− @C60 complex (∆∆E(DEF)(X = Li+ ) = 8.0, ∆∆E(DEF)(X = F− ) = −14.8 kcal/mol). However, the interaction energy increases when the Li+ ion is encapsulated (∆∆E(INT)A(BC)(X = Li+ ) = −35.9, ∆∆E(INT)A(BC)(X = F− ) = +33.8 kcal/mol), which shows that the encapsulated Li+ ion induces a stronger interaction between the metal fragment and the X@C60 fragment, so {η2 -(Li+ @C60)}PtL2

How Important is Metal-Carbon Back-Bonding for the Stability of Fullerene-Transition Metal Complexes?... http://dx.doi.org/10.5772/intechopen.70068 117

a Energies are given in kcal/mol.

**3.2. Basic energy decomposition analysis (basic EDA)**

) at the M06/LANL2DZ level of theory.

) − E(X@Cn) − E(ML2

**System Geometrical parameters NPA atomic charge**

116 Fullerenes and Relative Materials - Properties and Applications

ion induces more distortion in fragments A and B of the Pt-Li+

complexes, basic EDA is performed on {η2

ML2

Pt-F−

∆∆E(INT)A(BC)(X = F−

∆E= E({η2

M = Pt; X = Li+

M = Pt; X = 0

M = Pt; X = F−

M = Pd; X = Li+

M = Pd; X = 0

M = Pd; X = F−

ML2

ML2

ML2

ML2

ML2

ML2

({η2

The same results are obtained for X = Li+

@C60 complex (∆∆E(DEF)(X = Li+

interaction energy increases when the Li+

In order to better understand the factors that govern the thermodynamic stability of {η2

for n = 70, 76, 84, 90 and 96 are presented elsewhere). The bonding energy (∆E) is defined as

ions is initially considered. Basic EDA shows that both the metal fragment and the empty C60 fragment are distorted during the formation of the metal-carbon bond (**Table 2**). The metal fragment undergoes greater distortion (∆E(DEF)A = 52.1 kcal/mol) than C60 (∆E(DEF)B = 13.9 kcal/mol).

or F−

stronger interaction between the metal fragment and the X@C60 fragment, so {η2

) = 8.0, ∆∆E(DEF)(X = F−

) = +33.8 kcal/mol), which shows that the encapsulated Li+

**M-C1 M-C2 X-C1 X-C2 M C1 C2 X**

X@C60 2.12 2.12 2.29 2.29 +0.48 −0.32 −0.32 +0.86

X@C60 2.12 2.12 – – +0.46 −0.27 −0.27 –

X@C60 2.13 2.13 3.18 3.18 +0.43 −0.23 −0.23 −0.93

X@C60 2.12 2.12 2.28 2.28 +0.44 −0.27 −0.27 +0.86

X@C60 2.13 2.13 – – +0.40 −0.21 −0.21 –

X@C60 2.15 2.15 3.15 3.15 +0.34 −0.16 −0.16 −0.93

**Table 1.** Selected geometrical parameters (bond distances in Å) and the NPA atomic charge for optimized complexes

) = −35.9,

@C60)}PtL2

ion induces a

complexes (the basic EDA results

@C60 complex than those of the

) = −14.8 kcal/mol). However, the

), using Eq. (1). The Pt-C60 complex without encapsulated

ion is encapsulated (∆∆E(INT)A(BC)(X = Li+

. It is found that the encapsulation of the Li+

b A and B respectively represent the metal fragment and C60 cage.

c The difference is relative to corresponding quantity at X = 0.

dThe reaction energy without zero-point energy (ZPE) correction for the product, relative to the corresponding reactants.

**Table 2.** Basic EDA for {η2 -(X@C60)}ML2 (M = Pt, Pd) at M06/LANL2DZa,b .

is more stable. The relative thermodynamic stability increases in the order: ∆E(X = F− ) < ∆E (X = 0) < ∆E(X = Li+ ), as shown in **Table 2**. In addition, |∆∆E(DEF)| is small and |∆∆E(INT)A(BC)| is large, so the latter must be responsible for an increase in thermodynamic stability. Similar results are obtained for Pd-C60 complexes, but, the distortion in fragments A or B is smaller than that for the Pt-C60 complex, as is the interaction energy, so the complex is less stable. For instance, ∆E(M = Pd, X = Li+ ) = −34.5 > ∆E(M = Pt, X = Li+ ) = −35.9 kcal/mol. When the cage size increases (n = 70, 76, 84, 90 and 96), the encapsulated Li+ ion still induces a stronger interaction between the metal fragment and the X@Cn fragment and |∆∆E(DEF)| is smaller than |∆∆E(INT)A(BC)|. Therefore, an increase in the cage size has no effect on the basic EDA results.

#### **3.3. Advanced energy decomposition analysis (advanced EDA)**

In an earlier study by the authors, structural parameters and spectral characteristics were used to estimate the strength of π back-bonding for {η2 -(X@C60)}ML2 (M = Pt, Pd; X = 0, Li+ , L = PPh3 ) complexes [10]. The changes in bond length (Δr/r0 ), bond angle (Δθav), vibrational frequency (Δν) and the chemical shift (Δδ) were used to describe the character of the π-complex s. In this study, the strength of the π back-bonding strength is estimated from an energetic viewpoint using an advanced EDA method. This analysis shows the effect of encapsulated ions, metal fragments and cage sizes on π back-bonding.

#### *3.3.1. The effect of encapsulated ions on π back-bonding*

In an earlier discussion (Section 3.2), it was proven that thermodynamic stabilities increase in the order: ∆E(X = F− ) < ∆E(X = 0) < ∆E(X = Li+ ), because the interaction energy (∆E(INT)) is increased. The interaction between the metal fragment and X@Cn is now studied using advanced EDA, which further decomposes the interaction energy into electrostatic interaction (∆Eelstat), repulsive Pauli interaction (∆EPauli) and orbital interaction (∆Eorb) terms. The orbital interactions are the most important of these three and only the most important pairwise contributions to ΔEorb are considered. The advanced EDA method is used for {η2 -(X@ Cn)}ML2 complexes, as shown in **Tables 3** and **4** (the results for n = 70, 76, 84, 90 and 96 are presented elsewhere). A plot of the deformation density and a qualitative drawing of the orbital interactions between the metal fragment and X@C60 are shown in **Figure 2**. In terms of the Pt-C60 complexes, **Table 3** shows that both the electrostatic interaction (∆Eelstat) and the orbital interaction (∆Eorb) stabilize the complexes because they are negative terms, but the percentage of ∆Eorb increases in the order: ∆Eorb(X = F− ) < ∆Eorb(X = 0) < ∆Eorb(X = Li+ ). Therefore, the enhanced orbital interaction must be responsible for the increase in the thermodynamic stability. **Table 3** also shows that ΔE1 contributes significantly to ΔEorb: 69.5% for X = F− , 75.2% for X = 0 and 76.4% for X = Li+ . The deformation densities show that these come from π backdonation from a filled *d* orbital of the metal to the π\* orbitals of C60 (charge flow is yellow to green at the top of **Figure 2c**). The large contributions of ΔE1 to ΔEorb are in agreement with the results of previous studies. The metal-carbon bonds are principally formed by π backdonation [8]. It is also seen that the order of ΔE1 is |ΔE1 (X = F− )| = 94.4 < |ΔE1 (X = 0)| = 118.6 < |ΔE1 (X = Li+ )| = 142.8 kcal/mol. Therefore, ΔE1 is increased when there is the encapsulated Li+ ion but decreased when there is a F− ion. The second contribution of ΔE2 to ΔEorb is comparatively small: 18.0% for X = F− , 13.1% for X = 0 and 10.0% for X = Li+ . This results from

σ-donation from a filled π orbital of C60 to the π\* orbital of the metal (middle of **Figure 2c**). The computational results show that π back-bonding is crucial to the thermodynamic stability of

, 0, Li+

Pd-C60 complexes appear to be similar to Pt-C60 complexes, but a comparison of the results in

ing value for a Pt-C60 complex, which demonstrates that the π back-bonding for a palladium

 complexes (n = 60, 70, 76, 84, 90 and 96). It is seen that there is no linear relationship and there is one obvious peak for each X at n = 84 [26]. This demonstrates the effect of a difference in size of the carbon clusters on π back-bonding for a metal center, but the correlation is not simply monotonic. Therefore, a larger (smaller) cage size does not necessarily imply that there is stronger (weaker) π back-bonding, which results in greater (lower) thermodynamic

ion increases π back-bonding but an encapsu-

) at the M06/TZ2P level of theory. The fragments

for a Pd-C60 complex is smaller than the correspond-

)| = 142.8 kcal/mol. This is consistent with the earlier results that were

**Pd and C60 L2**

How Important is Metal-Carbon Back-Bonding for the Stability of Fullerene-Transition Metal Complexes?...

**Pd and Li+**

http://dx.doi.org/10.5772/intechopen.70068

**@C60**

119

values versus cage sizes that are calculated for {η2

(M = Pd, X = Li+

)| = 121.0

Pt-C60 complexes and that an encapsulated Li+

Optimized structures at the M06/LANL2DZ level of theory.

*3.3.2. The effect of metal fragments on π back-bonding*

center is weaker than that for a platinum center. For example, |ΔE1

obtained using structural parameters and spectral characteristics [10].

has the opposite effect.

**Table 4.** The advanced EDA results for {η2

**Tables 3** and **4** shows that the value ΔE1

*3.3.3. The effect of cage sizes on π back-bonding*

(M = Pt, X = Li+

**Figure 3** shows a plot of the ΔE1

lated F−

are PdL2

**Fragments L2**

ΔEelstat

ΔEorb

ΔE1

ΔE2

ΔE3

ΔE4

ΔE5

ΔErest

a

b

c

**Pd and F−**

**@C60 L2**

b −142.8 (61.0%) −137.0 (54.8%) −138.9 (49.6%)

b −91.2 (39.0%) −113.0 (45.2%) −141.3 (50.4%)

c −71.9 (78.8%) −96.9 (85.8%) −121.0 (85.6%)

c −15.5 (17.0%) −10.8 (9.6%) −9.4 (6.7%)

c −4.9 (5.4%) −4.2 (3.7%) −4.6 (3.3%)

c −2.1 (2.3%) −2.3 (2.0%) −3.6 (2.5%)

c −2.0 (2.2%) −2.9 (2.6%) −3.6 (2.5%)

c −2.7 (3.0%) −3.6 (3.2%) −6.5 (4.6%)

The values in parentheses give the percentage contribution to the total orbital interactions, ΔEorb.

and X@C60 in the singlet (S) electronic state. All energy values are in kcal/mol.

The values in parentheses give the percentage contribution to the total attractive interactions, ΔEelstat + ΔEorb.

a (X = F−

ΔEint −45.1 −73.5 −103.9 ΔEPauli 188.9 176.5 176.4

< |ΔE1

PtL2

stability.

a Optimized structures at the M06/LANL2DZ level of theory.

b The values in parentheses give the percentage contribution to the total attractive interactions, ΔEelstat + ΔEorb. c The values in parentheses give the percentage contribution to the total orbital interactions, ΔEorb.

**Table 3.** The advanced EDA results for {η2 -(X@C60)}PtL2 a (X = F− , 0, Li+ ) at the M06/TZ2P level of theory. The fragments are PtL2 and X@C60 in a singlet (S) electronic state. All energy values are in kcal/mol.

How Important is Metal-Carbon Back-Bonding for the Stability of Fullerene-Transition Metal Complexes?... http://dx.doi.org/10.5772/intechopen.70068 119

a Optimized structures at the M06/LANL2DZ level of theory.

is increased. The interaction between the metal fragment and X@Cn is now studied using advanced EDA, which further decomposes the interaction energy into electrostatic interaction (∆Eelstat), repulsive Pauli interaction (∆EPauli) and orbital interaction (∆Eorb) terms. The orbital interactions are the most important of these three and only the most important pairwise contributions to ΔEorb are considered. The advanced EDA method is used for {η2

presented elsewhere). A plot of the deformation density and a qualitative drawing of the orbital interactions between the metal fragment and X@C60 are shown in **Figure 2**. In terms of the Pt-C60 complexes, **Table 3** shows that both the electrostatic interaction (∆Eelstat) and the orbital interaction (∆Eorb) stabilize the complexes because they are negative terms, but the

the enhanced orbital interaction must be responsible for the increase in the thermodynamic

donation from a filled *d* orbital of the metal to the π\* orbitals of C60 (charge flow is yellow to

the results of previous studies. The metal-carbon bonds are principally formed by π back-

is |ΔE1

, 13.1% for X = 0 and 10.0% for X = Li+

percentage of ∆Eorb increases in the order: ∆Eorb(X = F−

118 Fullerenes and Relative Materials - Properties and Applications

donation [8]. It is also seen that the order of ΔE1

**Pt and F−**

Optimized structures at the M06/LANL2DZ level of theory.

**Table 3.** The advanced EDA results for {η2

ion but decreased when there is a F−

paratively small: 18.0% for X = F−

green at the top of **Figure 2c**). The large contributions of ΔE1

)| = 142.8 kcal/mol. Therefore, ΔE1

**@C60 L2**

ΔEint −67.7 −100.8 −133.0 ΔEPauli 257.2 235.8 227.6

b −189.0 (58.2%) −178.8 (53.1%) −173.8 (48.2%)

b −135.8 (41.8%) −157.7 (46.9%) −186.9 (51.8%)

c −94.4 (69.5%) −118.6 (75.2%) −142.8 (76.4%)

c −24.4 (18.0%) −20.6 (13.1%) −18.6 (10.0%)

c −7.1 (5.2%) −5.7 (3.6%) −6.5 (3.5%)

c −3.0 (2.2%) −3.8 (2.4%) −5.5 (3.0%)

c −3.6 (2.7%) −4.2 (2.7%) −5.0 (2.7%)

c – – −2.1 (1.1%)

c −5.8 (4.3%) −6.3 (4.0%) −6.8 (3.6%)

The values in parentheses give the percentage contribution to the total attractive interactions, ΔEelstat + ΔEorb.

a (X = F−

, 0, Li+

The values in parentheses give the percentage contribution to the total orbital interactions, ΔEorb.

and X@C60 in a singlet (S) electronic state. All energy values are in kcal/mol.

stability. **Table 3** also shows that ΔE1

for X = 0 and 76.4% for X = Li+

complexes, as shown in **Tables 3** and **4** (the results for n = 70, 76, 84, 90 and 96 are

) < ∆Eorb(X = 0) < ∆Eorb(X = Li+

)| = 94.4 < |ΔE1

is increased when there is the encapsulated

**Pt and Li+**

) at the M06/TZ2P level of theory. The fragments

contributes significantly to ΔEorb: 69.5% for X = F−

. The deformation densities show that these come from π back-

(X = F−

ion. The second contribution of ΔE2

**Pt and C60 L2**

Cn)}ML2

< |ΔE1

Li+

ΔEelstat

ΔEorb

ΔE1

ΔE2

ΔE3

ΔE4

ΔE5

ΔE6

ΔErest

a

b

c

are PtL2

(X = Li+

**Fragments L2**

). Therefore,

(X = 0)| = 118.6

to ΔEorb is com-

. This results from

**@C60**

to ΔEorb are in agreement with

, 75.2%

b The values in parentheses give the percentage contribution to the total attractive interactions, ΔEelstat + ΔEorb. c The values in parentheses give the percentage contribution to the total orbital interactions, ΔEorb.

**Table 4.** The advanced EDA results for {η2 -(X@C60)}PdL2 a (X = F− , 0, Li+ ) at the M06/TZ2P level of theory. The fragments are PdL2 and X@C60 in the singlet (S) electronic state. All energy values are in kcal/mol.

σ-donation from a filled π orbital of C60 to the π\* orbital of the metal (middle of **Figure 2c**). The computational results show that π back-bonding is crucial to the thermodynamic stability of Pt-C60 complexes and that an encapsulated Li+ ion increases π back-bonding but an encapsulated F− has the opposite effect.

#### *3.3.2. The effect of metal fragments on π back-bonding*

Pd-C60 complexes appear to be similar to Pt-C60 complexes, but a comparison of the results in **Tables 3** and **4** shows that the value ΔE1 for a Pd-C60 complex is smaller than the corresponding value for a Pt-C60 complex, which demonstrates that the π back-bonding for a palladium center is weaker than that for a platinum center. For example, |ΔE1 (M = Pd, X = Li+ )| = 121.0 < |ΔE1 (M = Pt, X = Li+ )| = 142.8 kcal/mol. This is consistent with the earlier results that were obtained using structural parameters and spectral characteristics [10].

#### *3.3.3. The effect of cage sizes on π back-bonding*

**Figure 3** shows a plot of the ΔE1 values versus cage sizes that are calculated for {η2 -(X@Cn)} PtL2 complexes (n = 60, 70, 76, 84, 90 and 96). It is seen that there is no linear relationship and there is one obvious peak for each X at n = 84 [26]. This demonstrates the effect of a difference in size of the carbon clusters on π back-bonding for a metal center, but the correlation is not simply monotonic. Therefore, a larger (smaller) cage size does not necessarily imply that there is stronger (weaker) π back-bonding, which results in greater (lower) thermodynamic stability.

**4. Conclusion**

stability of {η2

F−

a F−

between cage size and π back-bonding.

ion has the opposite effect.

**Figure 3.** The correlation between ΔE1

blue, red and black lines indicate the ΔE1

**Acknowledgements**

but the complex becomes unstable if there is a F−

This computational study uses density functional theory to determine the thermodynamic

How Important is Metal-Carbon Back-Bonding for the Stability of Fullerene-Transition Metal Complexes?...

, respectively.

The advanced EDA results show that π back-bonding is crucial to thermodynamic stability

The authors would like to thank the National Center for High-Performance Computing in Taiwan for the donation of generous amounts of computing time. The authors are also grate-

ful for financial support from the Ministry of Science and Technology of Taiwan.

 ion has the opposite effect. These computations also show that a platinum center results in stronger π back-bonding than a palladium center and that there is no linear relationship

, 0, Li+

and n = 60, 70, 76, 84, 90 and 96).

(n = 60, 70, 76, 84, 90 and 96) complexes. The

http://dx.doi.org/10.5772/intechopen.70068

121

ion is encapsulated within Cn

ion but the presence of

ion but

ion. Basic EDA shows that there is an increase

complexes (M = Pt, Pd; X = F−

and cage sizes for {η2

, 0 and F−

for X = Li+

in the interaction between the metal fragment and Cn if there is an encapsulated Li+

The calculations show the reaction is more stable when the Li+

and that thermodynamic stability is increased by the presence of a Li+

**Figure 2.** (a) A qualitative drawing of the orbital interactions between the metal fragment and Li+ @C60; (b) the shape of the most important interacting occupied and vacant orbitals of the metal fragments and Li+ @C60; (c) a plot of the deformation densities, Δρ, for the pairwise orbital interactions between the two fragment in their closed-shell state, the associated interaction energies, ΔEorb (in kcal/mol), and the eigenvalues ν. The eigenvalues, ν, indicate the size of the charge flow. The direction of the charge flow is from yellow to the green.

How Important is Metal-Carbon Back-Bonding for the Stability of Fullerene-Transition Metal Complexes?... http://dx.doi.org/10.5772/intechopen.70068 121

**Figure 3.** The correlation between ΔE1 and cage sizes for {η2 -(X@Cn)}PtL2 (n = 60, 70, 76, 84, 90 and 96) complexes. The blue, red and black lines indicate the ΔE1 for X = Li+ , 0 and F− , respectively.