Rachel O Mariani

Marc W Cadotte

Marney E Isaac

Denis Vile

Cyrille Violle

Adam R Martin
email: adam.martin@utoronto.ca

National-scale changes in crop diversity through the Anthropocene

Expansion of crops beyond their centres of domestication is a defining feature of the Anthropocene Epoch. This process has fundamentally altered the diversity of croplands, with likely consequences for the ecological functioning and socio-economic stability of agriculture under environmental change. While changes in crop diversity through the Anthropocene have been quantified at large spatial scales, the patterns, drivers, and consequences of change in crop diversity and biogeography at nationalscales remains less explored. We use production data on 339 crops, grown in over 150 countries from 1961 to 2017, to quantify changes in country-level crop richness and evenness. Virtually all countries globally have experienced significant increases in crop richness since 1961, with the early 1980s marking a clear onset of a ~ 9-year period of increase in crop richness in countries worldwide. While these changes have increased the similarity of diversity of croplands among countries, only half of countries experienced increases in crop evenness through time. Ubiquitous increases in crop richness within nearly all countries between 1980 and 2000 are a unique biogeographical feature of the Anthropocene. At the same time, we detected opposing changes in crop evenness, and only modest signatures of increased homogenization of croplands among countries. Therefore context-dependent and, at least, national-scale assessments are needed to understand and predict how changes in crop diversity influence agricultural resistance and resilience to environmental change.

The Anthropocene Epoch defines the current epoch of geologic time, where human activity is the dominant force shaping Earth's abiotic and biotic environmental processes, systems, and cycles [START_REF] Crutzen | Earth System Science in the Anthropocene[END_REF][START_REF] Lewis | Defining the Anthropocene[END_REF][START_REF] Malhi | The concept of the Anthropocene[END_REF] . This novel geological epoch is defined in part by human-caused changes to Earth's biogeography e.g. 4 , including the introduction of non-native or domesticated species e.g. [START_REF] Capinha | The dispersal of alien species redefines biogeography in the Anthropocene[END_REF] , and climate-and land-use change-induced shifts in species distributions e.g. [START_REF] Parmesan | A globally coherent fingerprint of climate change impacts across natural systems[END_REF] . However, arguably the most pronounced human-mediated changes in species biogeography defining the Anthropocene, albeit still largely understudied [START_REF] Zimmerer | The biodiversity of food and agriculture (Agrobiodiversity) in the anthropocene: Research advances and conceptual framework[END_REF] , is the deliberate spread of crops outside of their centres of domestication into other parts of the world [START_REF] Martin | Regional and global shifts in crop diversity through the Anthropocene[END_REF][START_REF] Crosby | The Columbian Exchange: Biological and Cultural Consequences of 1492[END_REF][START_REF] Mann | How the ecological collision of Europe and the Americas gave rise to the modern world[END_REF] . Large-scale anthropogenic influences on crop biogeography first emerged during the Columbian Exchange, which is a period in the early sixteenth century defined by massive exchanges of crop plants and animals between West Africa, Europe, and the Americas, during the time of New World colonization by Europeans [START_REF] Mann | How the ecological collision of Europe and the Americas gave rise to the modern world[END_REF] . During this time there was major movement of crops such as corn and beans from the Americas to Europe, while other crops such as wheat and barley were transported the opposite direction [START_REF] Crosby | The Columbian Exchange: Biological and Cultural Consequences of 1492[END_REF][START_REF] Mann | How the ecological collision of Europe and the Americas gave rise to the modern world[END_REF] .

While research defining the scientific basis of the Anthropocene have focused on the Columbian Exchange 2 , analyses of contemporary global food production trends have demonstrated that there have been more recent (i.e., post-1950) and widespread movement of new crops-i.e., those previously not present in large, industrial agricultural lands-into nearly all regions of the world [START_REF] Martin | Regional and global shifts in crop diversity through the Anthropocene[END_REF] . Specifically since the early 1960s, the range of crops cultivated in agricultural lands across virtually all continents and regions, has changed in remarkably similar patterns: (1) in the 1960s continents experienced a period of little change in the number of new crops being cultivated in large (i.e., industrial) agricultural lands; this was followed by (2) a period of rapid increase in the number of crop groups being cultivated, beginning in the late 1970s and continuing through to the 1980s; and finally, (3) there exists a period of little change or levelling off in the number of crops being cultivated beginning in the 1990s and persisting through to present day [START_REF] Martin | Regional and global shifts in crop diversity through the Anthropocene[END_REF] . These changes have drastically shaped not only the diversity of agricultural lands worldwide, but have greatly influenced global food supplies and security, diets, and agricultural economics [START_REF] Aguiar | Global changes in crop diversity: Trade rather than production enriches supply[END_REF][START_REF] Campi | Specialization in food production affects global food security and food systems sustainability[END_REF][START_REF] Khoury | Origins of food crops connect countries worldwide[END_REF][START_REF] Khoury | Increasing homogeneity in global food supplies and the implications for food security[END_REF][START_REF] Nelson | Commercial plant production and consumption still follow the latitudinal gradient in species diversity despite economic globalization[END_REF][START_REF] Renard | National food production stabilized by crop diversity[END_REF][START_REF] Mahaut | Complementary mechanisms stabilize national food production[END_REF] .

Coinciding with these human-mediated increases in the number of crops being cultivated in agricultural lands, this time period (particularly the 1980s) is also marked by an overall increase in the similarity of crop composition globally [START_REF] Martin | Regional and global shifts in crop diversity through the Anthropocene[END_REF] . As an example of this homogenization of croplands at supra-national scales, wheat, maize, soy, and rice, now dominate over 50% of the world's agricultural lands 18 . These changes have contributed to major shifts in the world's agricultural landscapes and economies, including the homogenization of global food supplies and diets [START_REF] Khoury | Increasing homogeneity in global food supplies and the implications for food security[END_REF] , greater interdependency in agricultural trade between countries to maintain food security and potentially increases susceptibility of agricultural lands to pests, diseases, and climate stressors [START_REF] Hessenauer | Evolution and adaptation of forest and crop pathogens in the Anthropocene[END_REF] .

Researchers have argued that the magnitude of these contemporary changes in crop biogeography, indicates that the 1980s are also a notable marker of the onset of the Anthropocene 8 . However, this aspect of a crop biogeography-based line of evidence supporting the Anthropocene hypothesis remains limited to research on changes in crop diversity at very large supranational regional, continental, or global spatial scales [START_REF] Martin | Regional and global shifts in crop diversity through the Anthropocene[END_REF] , or a small number of national-scale analyses (e.g., within the United States [START_REF] Aguilar | Crop species diversity changes in the United States: 1978-2012[END_REF] ). This is despite reason to expect that the timing and rates of change in crop diversity is likely to differ widely among countries. For example previous work has shown that regional-scale agricultural policy initiatives including the Caribbean Basin Initiative and the North American Free Trade Agreement, were a primary catalyst for changes in crop diversity in the Caribbean and North American, respectively [START_REF] Martin | Regional and global shifts in crop diversity through the Anthropocene[END_REF][START_REF] Isakson | Maize diversity and the political economy of agrarian restructuring in Guatemala[END_REF] . Yet such regional initiatives are not necessarily common in many agricultural economies. For instance, throughout regions of Africa, country-specific structural readjustment policies and colonial histories, but not necessarily region-scale policy initiatives, have played a key role in determining crop composition on agricultural lands [START_REF] Bonneuil | Development as experiment: Science and state building in late colonial and postcolonial Africa, 1930-1970[END_REF] . Similarly, in certain countries Western philanthropic organizations have played a significant role in transferring crop technologies including disease resistant crops and pressuring shifts to cash crops, which have had unintended consequences for national agrobiodiversity [START_REF] Bonneuil | Development as experiment: Science and state building in late colonial and postcolonial Africa, 1930-1970[END_REF][START_REF] Holt-Gimenez | Food First Policy Brief No. 12: Ten Reasons Why the Rockefeller and the Bill and Melinda Gates Foundations' Alliance for Another Green Revolution Will Not Solve the Problems of Poverty and Hunger in Sub-Saharan Africa[END_REF] . This and a multitude of similar examples (e.g., [START_REF] Isakson | Maize diversity and the political economy of agrarian restructuring in Guatemala[END_REF][START_REF] Arnold | technology, and colonialism in the 20th century[END_REF][START_REF] Meertens | Agricultural performance in Tanzania under structural adjustment programs: Is it really so positive?[END_REF][START_REF] Chen | Diversification and intensification of agricultural adaptation from global to local scales[END_REF] ) indicate that country-scale analyses are needed for a nuanced understanding of how crop diversity and composition has changed through the Anthropocene.

Supranational assessments also preclude nuanced analyses of the correlates or consequences of crop diversity change. Specifically, much of the structural adjustment policies that emphasized enhanced balance of trade and exports were imposed on the agricultural economies of developing nations throughout the 1980s e.g. [START_REF] Isakson | Maize diversity and the political economy of agrarian restructuring in Guatemala[END_REF][START_REF] Meertens | Agricultural performance in Tanzania under structural adjustment programs: Is it really so positive?[END_REF] . Since these policies tended to focus on cultivation of new crops for international export markets (e.g., the introduction of cocoa in India and oil palm in Peru in the early and mid-1970s), one might expect that changes in agricultural diversity profiles (i.e., the timing, duration, and rate of change in crops being cultivated) vary systematically with socio-economic development indices. More specifically, one may hypothesize that a certain socio-economic group (e.g., developing nations with lower Human Development Index (HDI) values) to have broadly shown similar timing and rates of change in crop diversity, compared to countries with higher HDI.

Finally, one might expect that the number of crops cultivated within a given country also varies systemically with latitude, across a latitudinal gradient in crop diversity [START_REF] Nelson | Commercial plant production and consumption still follow the latitudinal gradient in species diversity despite economic globalization[END_REF] . This pattern is likely to emerge as a function of multiple factors including (1) seasonality and limited growing conditions towards the poles [START_REF] Pontarp | The latitudinal diversity gradient: Novel understanding through mechanistic eco-evolutionary models[END_REF] , and (2) the centres of crop domestication being disproportionately situated in tropical and sub-tropical regions [START_REF] Harlan | Agricultural origins: Centers and noncenters[END_REF][START_REF] Meyer | Patterns and processes in crop domestication: An historical review and quantitative analysis of 203 global food crops[END_REF] . Moreover, if a latitudinal gradient in the number of crops cultivated in agricultural lands does exist, these patterns may have been drastically altered by the incorporation of new crops into agricultural lands that began in the 1980s [START_REF] Martin | Regional and global shifts in crop diversity through the Anthropocene[END_REF] .

National-scale analyses of crop cultivation are needed to test these questions surrounding crop diversity change, and its role in defining the Anthropocene. To address this, we used crop production data collected by the Food and Agricultural Organization (FAO) 18 to execute national-scale analyses of changes in crop diversity. This analysis evaluated production data for 339 crop species groups, grown in 201 United Nations-recognized countries, over the past 56 years (from 1961 to 2017), resulting in a dataset of over ~ 2.36 million data points. Our analysis was designed to address the following: (1) How do patterns of change in the diversity of crops grown on agricultural lands, vary across countries in recent decades? (2) Have changes in country-scale crop diversity resulted in a detectable signal of "homogenization" across the world's agricultural lands? (3) Are patterns of change in crop diversity across countries, systematically related to country-scale socio-economic indicators? (4) Does the number of crops cultivated across countries follow a latitudinal diversity gradient? And if so, (5) has the introduction of new crops into different countries altered this gradient?

Results

Changes in crop richness within countries through time. Our results are based on trends and three different "Indicator" values inferred from a piecewise modelling framework, which is described in detail in the Methods and presented conceptually in Fig. 1. Across all countries we detected significant changes, largely increases, in crop commodity group richness over time (r 2 = 0.356-0.998 across 165 countries, p < 0.01 in all cases; Fig. 2, Table S1). Across 165 countries for which piecewise models converged, the average initial onset of crop group richness changes (Indicator 1) occurred in 1983 (± 9.2 years s.d.; Fig. 2A,B). The onset of changes in crop group richness occurred at the earliest in 1962 (in India), with the most recent onset of crop group richness change beginning in 2007 (in Serbia; Fig. 2A,B). Richness began to change at or after 1980 in 113 countries, of which 71% (or 80 countries) show changes in crop richness beginning from 1980 to 1989 (Fig. 2A; Table S1).

The period of crop richness change (Indicator 2) lasted approximately 9 years on average (± 8.9 years m.a.d.; Fig. 2C,D). Nineteen countries had a period of richness change lasting only one year, while nine countries experienced a more prolonged period of crop richness change lasting ≥ 30 years (Fig. 2C,D). National changes in richness began to saturate (i.e., ψ2 in Eq. ( 2), or breakpoint 2 in Fig. 1) on average in 1995 (± 8.3 years s.d.), with richness increases levelling off in 117 countries prior to the turn of the century at or before 1999 (Fig. 2C,D).

Over the period of change (i.e., between ψ1 and ψ2, Fig. 1), crop richness increased in 151 of 165 countries evaluated. Across all countries, richness increased (Indicator 3) on average by 0.8 species per year (± 0.9 m.a.d.; Fig. 2E,F). Only 13 countries reported decreases in richness per year, with these countries showing declines of 0.4 ± 0.6 crop groups per year on average (Fig. 2E,F, Table S1). In all of our datasets, Indicators 1-3 associated with changes in S were unrelated to the total area of a country (linear regression r 2 = 0.0-0.004, p ≥ 0.888), or the area of a country under crop cultivation (linear regression r 2 = 0.0-0.01, p ≥ 0.208).

Across the 164 countries for which data was available in both 1961 and 2017, crop richness varied significantly as a function of latitude in similar patterns in both years (Table 1, Fig. 3). Latitude and a 2nd order polynomial term explained 10.5% and 12.9% of the variation in crop group richness in 1961 and 2017, respectively (model p ≤ 0.001 in both cases), with the richness-latitude relationship being similar in both years (Table 1). Specifically, regression models indicated that crop group richness increased from equatorial regions towards mid-latitude countries, with modeled peak crop group richness occurring at ~ 37° latitude in both 1961 and 2017; richness then declined at higher latitudes, denoted by statistically significant (p ≤ 0.002) negative 2nd order polynomial terms (Fig. 3). While these trends were similar between years, this analysis did reflect the increased crop group richness that occurred between 1961 and 2017, in 152 of 164 of the countries included in this analysis.

Changes in crop evenness and composition across countries through time.

Piecewise models evaluating changes in crop group evenness (J′) through time converged for 185 countries; in these countries year explained on average 80.9% of the variation in J′ (model p < 0.01 in all cases, r 2 range = 0.198-0.992; Table S2). In these countries changes in J′ began on average in 1981 (± 11.9 years m.a.d.), although compared to changes in group richness, periods of change in evenness were more prolonged lasting on average for 15 years (± 13.3 years m.a.d.) (Fig. 4A-D). Unlike analyses of crop richness, changes in evenness were less systematic, such that through the period of change evenness declined in 97 countries, increased in 88 countries, and average changes in J′ centred on zero (mean = -0.004 year -1 ± 0.04 s.d.; Fig. 4E,F). However, similar to patterns of change in richness, Indicators 1-3 associated with changes in evenness were also independent of total area of a country Values for all three indicators for each country were derived from piecewise linear models predicting S as a function of year (see Fig. 1 for example). Countries coloured gray in the maps were those where either data was not available or the piecewise models failed to converge (denoted in Table S1). Histograms and associated descriptive statistics for each indicator are also presented, with means (± s.d.) or medians (± m.a.d.) denoted visually by the points and error bars below the histograms. All piecewise model parameters for each country are presented in Table S1. Maps were generated using the mapCountryData function in the rworldmap R package [START_REF] South | rworldmap: A new R package for mapping global data[END_REF] .

(linear regression r 2 = 0.001-0.004, p ≥ 0.404), or the area of a country under crop cultivation (linear regression r 2 = 0.0-0.013, p ≥ 0.12).

Multivariate analysis detected a significant influence of country, year, and a year-by-country interaction term on crop composition (Adonis p < 0.01 in all cases; Table 2). Of these variables, country differences were most pronounced with country identity explaining 89.5% of the variability in crop composition. Year explained ~ 1% of the variability in crop composition, while a country-by-year interaction term explained an additional 6.3% of the variation in crop composition (Table 2). Based on the non-metric multidimensional scaling (NMDS) analysis, there was trend of increasing similarity in crop composition among countries over time. This is illustrated by increasingly smaller 95% confidence bands surrounding the data points along NMDS axes 1 and 2 from 1961 through to 1983 (i.e., approximately the average year in which changes in crop group richness and evenness commence); multivariate space encapsulated by the 95% confidence band is then further reduced through 2017 (Figure S1). Table 1. Parameters and diagnostics for linear regression models (including a 2nd-order polynomial term) predicting crop group richness as a function of latitude in 1961 and 2017. Only the 164 countries with data from both years were included in these analyses. Parameter estimate standard errors shown in parentheses. Also shown are the latitudes at which crop group richness peaked, according to these models, which are shown visually in Fig. 3. Only countries with data from both years are included in this analysis. Trend lines represent results from linear models (including a 2nd-order polynomial term) that predict crop group richness as a function of latitude (and latitude 2 ) in 1961 (dashed line, filled circles) and 2017 (solid line, open circles). Complete diagnostics for both models are presented in Table 1.

Year

Socioeconomics correlates of changes in crop group diversity and evenness. Patterns of change

in crop richness or evenness were not systematically related to country socio-economic status, with HDI values predicting only 1.3% of Indicators 1-3 for both S and J′ on average (Table 3). The only exception to this is that Values for all three indicators for each country were derived from piecewise linear models predicting J′ as a function of year, where harvested area (in ha) was used to approximate group abundance. Countries coloured gray in the maps were those where either data was not available or the piecewise models failed to converge (denoted in Table S2). Histograms and associated descriptive statistics for each indicator are also presented, with means (± s.d.) or medians (± m.a.d.) denoted visually by the points and error bars below the histograms. All piecewise model parameters for each country are presented in Table S2. Maps were generated using the mapCountryData function in the rworldmap R package [START_REF] South | rworldmap: A new R package for mapping global data[END_REF] . the rate of change in S was significantly negatively related to HDI: countries with lower HDI scores expressed greater increases in S (i.e., Indicator 3 values) during their period of crop diversity change vs. countries with higher HDI (Table 3). This analysis did point to stronger explanatory power of spatial location in determining patterns of change in crop S and J′: both continent and region identity explained an average of 9.5% and 2.0% of the variation in Indicators 1-3, respectively (Table 3). However, ultimately country-by-country variation in crop group diversity change was largely idiosyncratic, with 87.0% of the variation in Indicators 1-3 for both S and J′ being unaccounted for by socio-economic and spatial factors included in our models.

Discussion

The vast majority of the world's countries have experienced significant increases in the number of crop groups being cultivated over recent decades; changes that contribute to a detectable increase in the similarity of crops being grown among countries, that has occurred since the 1960s. Our findings align with those from previous studies [START_REF] Martin | Regional and global shifts in crop diversity through the Anthropocene[END_REF] that detected similar changes in crop taxonomic and phylogenetic diversity at supra-national regional, continental, and global scales: (1) a period of little change through the 1960s and 1970s, followed by (2) the onset of increases in crop commodity group richness commencing the 1980s and extending through the end of the 1990s, followed by (3) a levelling off of commodity group richness beyond the 2000s. While the concordance of these two studies should not be surprising, our findings contribute a new and more nuanced understanding of the remarkably similar patterns in crop group richness that have occurred across virtually all nations in recent decades; the consistency of which supports the idea that changes in crop diversity and biogeography occurring since the Columbian Exchange, are a ubiquitous feature of the Anthropocene.

However, these near-universal patterns in crop richness increases over recent decades (observed in all but 13 countries evaluated here) have had unequal and more variable effects on the evenness of crops being cultivated. The average duration of change in evenness (15 years) was nearly twice as long vs. the average duration of changes in richness (9 years). Previous work has noted that the majority of introduced crops into national agricultural production portfolios, largely owes to cultivation of crops beyond their country or region of origin [START_REF] Khoury | Origins of food crops connect countries worldwide[END_REF] . So while introduction of novel crop groups has been rather succinct-corresponding to periods of rapid change in structural adjustment through the 1980s and 1990s [START_REF] Isakson | Maize diversity and the political economy of agrarian restructuring in Guatemala[END_REF][START_REF] Bonneuil | Development as experiment: Science and state building in late colonial and postcolonial Africa, 1930-1970[END_REF] -expansion of these crops across more cultivated lands draws out over longer periods. This discrepancy likely reflects the lag between new crop introductions, compared to longer-term process associated with agricultural economic adjustments for these crops including expansion of export markets or crop-specific subsidies [START_REF] Nelson | Commercial plant production and consumption still follow the latitudinal gradient in species diversity despite economic globalization[END_REF] , and to a much lesser extent, expansion of domestic consumption markets for new crops [START_REF] Aguiar | Global changes in crop diversity: Trade rather than production enriches supply[END_REF] .

The specific model parameters for country-specific trends may be sensitive to data quality [START_REF] Mahaut | Complementary mechanisms stabilize national food production[END_REF] . However, here we observed a clear signal of both decreasing and increasing trends in crop evenness across 52% and 48% of the world's countries, respectively. This approximately even proportion of increase and decrease in evenness detected among all countries here, would explain why analyses of global scale crop production trends reported no change in evenness between 1961 and 2013 [START_REF] Aguiar | Global changes in crop diversity: Trade rather than production enriches supply[END_REF] . Yet similar to the growing number of analyses of production and consumption at multiple scales [START_REF] Martin | Regional and global shifts in crop diversity through the Anthropocene[END_REF][START_REF] Aguiar | Global changes in crop diversity: Trade rather than production enriches supply[END_REF][START_REF] Khoury | Increasing homogeneity in global food supplies and the implications for food security[END_REF][START_REF] Nelson | Commercial plant production and consumption still follow the latitudinal gradient in species diversity despite economic globalization[END_REF] , our multivariate analysis does indicate that over recent decades, countries have expressed a statistically significant increase in the similarity of crops in their agricultural lands. The implications of these shifts remain speculative, and based on our analysis, are clearly scale-dependent. Specifically, increases in the similarity in crop composition across regions and continents could indicate growing susceptibility of agriculture to pest and pathogen outbreaks, and perhaps climate change effects like regional temperature increases or precipitation declines [START_REF] Martin | Regional and global shifts in crop diversity through the Anthropocene[END_REF] . Alternatively, a wider geographic spread of crop groups could well represent a means of buffering production from localized disruptions including local climatic change or weather events, pest or pathogens, or civil unrest [START_REF] Aguiar | Global changes in crop diversity: Trade rather than production enriches supply[END_REF] . While both are plausible, generalizing either hypothesis to predict the impacts of homogenization across all agricultural lands globally is inconsistent with our countryspecific analyses here. Since 1961, the evenness of crops in agricultural lands is both increasing and decreasing in approximately the same proportion of countries and agricultural area globally. Therefore, while previous studies including our own have speculated on how changes in crop diversity will influence global agricultural production and sustainability [START_REF] Martin | Regional and global shifts in crop diversity through the Anthropocene[END_REF][START_REF] Aguiar | Global changes in crop diversity: Trade rather than production enriches supply[END_REF] , our results suggest the largest spatial scale at which potential impacts of stressors on crop production-though not food consumption or food security per se-can be robustly predicted is on a per-country basis.

Knowledge of the high context-dependency of agricultural adaptation, management, and crop selection likely indicates even smaller spatial scales (i.e., communities, households) are needed in order to fully predict susceptibility of production in the future. Indeed, an important caveat for our analysis is that even national scales likely do not comprehensively indicate how agricultural functional diversity has changed in the past 60 years. Our work here focuses on large-scale industrial agriculture, with the contributions of small-scale agriculture to overall crop diversity, particularly in terms of locally adapted varieties or landraces, likely underestimated [START_REF] Ricciardi | Higher yields and more biodiversity on smaller farms[END_REF] . Spatially explicit records of crop genetic diversity including locally adapted or cultivated crop phenotypes and their wild relatives are becoming more widely available, yet such a comprehensive assessment of crop genetic diversity remains prohibitive; indeed, while instructive, even analyses of phylogenetic crop diversity [START_REF] Martin | Regional and global shifts in crop diversity through the Anthropocene[END_REF][START_REF] Milla | Crop origins and phylo food: A database and a phylogenetic tree to stimulate comparative analyses on the origins of food crops[END_REF][START_REF] Milla | Phylogenetic patterns and phenotypic profiles of the species of plants and mammals farmed for food[END_REF] are limited to the crop species or sub-species levels. Clearly more work is needed to better integrate the methodological frameworks and concepts used here and in related studies [START_REF] Aguiar | Global changes in crop diversity: Trade rather than production enriches supply[END_REF][START_REF] Khoury | Increasing homogeneity in global food supplies and the implications for food security[END_REF] , with finer-scale datasets on crop genetic diversity.

Determination of potential driving factors behind diversification trends. The factors underpin-

ning national-scale increases in crop richness or evenness include a nuanced mix of environmental, socio-economic, and cultural factors. Here we hypothesized that patterns of change in crop group richness would correlate with development status (quantified as the Human Development Index). However, HDI was clearly a poor predictor of the rate, duration, and timing of changes in crop group richness and evenness. While national-scale agriculture portfolios did not systematically change as a function of HDI, spatial location (i.e., region and continent) did explain ~ 10-12% of the variation in patterns of change. From a strictly socio-economics perspective, this finding could point to regional-scale agricultural policy in driving similar changes in crop group richness, on a region-by-region basis [START_REF] Isakson | Maize diversity and the political economy of agrarian restructuring in Guatemala[END_REF][START_REF] Meertens | Agricultural performance in Tanzania under structural adjustment programs: Is it really so positive?[END_REF] . Alternatively, it could also reflect regional-scale similarities in non-governmental organization (NGO) interventions in the agricultural sector [START_REF] Chen | Diversification and intensification of agricultural adaptation from global to local scales[END_REF] . Indeed, previous studies have indicated that governmental and NGO intervention has led to a ~ 17% increase in global crop diversity on average [START_REF] Chen | Diversification and intensification of agricultural adaptation from global to local scales[END_REF] . Additionally, and perhaps surprisingly, this same study [START_REF] Chen | Diversification and intensification of agricultural adaptation from global to local scales[END_REF] found that while climate is an important driver of on-farm crop diversity change, its influence is secondary compared to market conditions.

Agriculture and the reshuffling of species through the Anthropocene. We detected a statistically significant and hump-shaped latitudinal gradient in crop group richness, with low-to mid-latitude regions (centred on ~ 37° latitude), expressing the largest number of cultivated crop groups in large-scale agricultural lands. Consistent with previous studies, this non-linear latitudinal trend most likely emerges due to: (1) countries at these latitudes having a range of climatic conditions (i.e., Köppen climate zones) that supports cultivation of a large diversity of crop functional types and year-round production; and (2) countries at these latitudes encapsulating many of the world's centres of crop domestication, and having received among the largest imports of crops during the Columbian Exchange [START_REF] Khoury | Origins of food crops connect countries worldwide[END_REF][START_REF] Meyer | Patterns and processes in crop domestication: An historical review and quantitative analysis of 203 global food crops[END_REF][START_REF] Milla | Phylogenetic patterns and phenotypic profiles of the species of plants and mammals farmed for food[END_REF] .

It is expected that the band of latitude supporting the highest number of crop groups could shift as a result of global environmental change drivers. Indeed, researchers have now long projected that crop diversity and richness within countries will change as climate change intensifies, with clusters of high crop diversity moving poleward through time [START_REF] Leemans | Modeling the potential change in yield and distribution of the earth's crops under a warmed climate[END_REF][START_REF] Zabel | Global agricultural land resources-a high resolution suitability evaluation and its perspectives until 2100 under climate change conditions[END_REF] . However, similar to categorical comparisons of crop richness produced in tropical vs. temperature countries [START_REF] Nelson | Commercial plant production and consumption still follow the latitudinal gradient in species diversity despite economic globalization[END_REF] , our analysis did not detect evidence of a disruption in the latitudinal trends of crop species richness in between 1961 and 2017. Instead, we find only a systematic increase in richness across all latitudes, indicating that the cultivation of new crop commodity groups in large-scale agricultural lands has not fundamentally altered the latitudes at which crop group richness is highest. Moreover, our analysis here (1) likely misses the role that high crop diversity of small holder (i.e. < 2-ha in size) farms [START_REF] Jarvis | A global perspective of the richness and evenness of traditional crop-variety diversity maintained by farming communities[END_REF] -which are disproportionately concentrated at lower latitudes [START_REF] Lowder | The number, size, and distribution of farms, smallholder farms, and family farms worldwide[END_REF] -play in driving latitudinal crop diversity gradients; and (2) does not address functional-or phylogenetic crop diversity, which may show different latitudinal patterns. Including these factors in additional analyses is a key step for further resolving our understanding of how crop biogeography is changing during the Anthropocene.

Methods

Data acquisition. Our analysis was based on open access crop production data from the United Nation's Food and Agricultural Organization (FAO) spanning from 1961 to 2017 18 . We extracted data on area harvested (in ha) for 339 FAO-defined crop groups being grown in all UN-recognized countries. Since our research centred on understanding, quantifying, and mapping changes in crop diversity in current agricultural lands, countries that cease to exist (e.g., Yugoslavia) were not included in our analysis, resulting in data for 201 countries (Table S1). Prior to analyses, we adjusted certain crop group listings following our previous analyses of global changes in crop diversity [START_REF] Martin | Regional and global shifts in crop diversity through the Anthropocene[END_REF] . Specifically, "Cottonlint" and "Cottonseed" were duplicated in our dataset and were therefore compiled as "Seedcotton", while "Palmkernels" were renamed as "Oilpalmfruit. " Additionally, "Fruitpomenes", "Fruitstonenes", and "Grainmixed" were removed from analysis since these crop groupings are not associated with any specific crop species in the FAO database 18 . Finally, "Mushroomsandtruffles" were removed since it relates to non-plant species, and "Coir" was removed because it is a plant by-product.

Changes in crop richness over time. All statistical analyses were performed using R version 3.3.3 statistical software (R Foundation for Statistical Computing, Vienna, Austria). The initial step in our analysis was to calculate both crop richness and evenness for each country, at each individual year, using the vegan R package [START_REF] Oksanen | Community Ecology Package. R package version 2[END_REF] . Based on these datasets, we then used the analytical framework developed by [START_REF] Martin | Regional and global shifts in crop diversity through the Anthropocene[END_REF] to evaluate how crop species richness and evenness have changed in each individual country across its entire data range.

Specifically, in their analysis Martin et al. [START_REF] Martin | Regional and global shifts in crop diversity through the Anthropocene[END_REF] found that piecewise linear regression models provided the strongest descriptions of crop species richness change over time, across 21 of 22 FAO-defined regions globally. We therefore followed this approach by fitting a piecewise linear regression model for each country individually, that predicts changes in species richness over time. Piecewise model fitting was a two-step process, whereby for each country we first fit a linear regression model of the form:

where a is the intercept and b represents the rate of change in crop group richness (S) through time. This linear model (Eq. 1) was then used as the basis of a piecewise linear regression model, which was fitted in order to estimate breakpoints in the relationship between S and year. Specifically, piecewise models were fit using the segmented function in the segmented R package [START_REF] Vito | Segmented: An R package to fit regression models with broken-line relationships. cran.r-project[END_REF] , and were of the form:

where a is as in Eq. ( 1), and b represents the slope of the S-year relationship prior to the first breakpoint (ψ1). Here, c represents the difference in the slope of the S-year relationship between the first and second piecewise model segments; the c parameter therefore applies only when the first conditional indicator function (denoted by "I") is true. Similarly, d represents the difference in slopes for the S-year relationship between the first, second, and third segments, which only applies when the second conditional indicator function is true. In sum, the slope of the relationship between S and year is equal to b prior to the ψ1, is equal to b + c between ψ1and ψ2, and is equal to b + c + d after ψ2. Piecewise models were fit with initial starting parameters of 1975 and 2000 for ψ1 and ψ2, respectively. The ψ1 and ψ2 parameters were tuned manually for 29 countries with a shortened data range, following visual inspection of data (see Tables S1 andS2).

Based on this piecewise regression model procedure, we then used parameters from Eq. (2) to determine three key indicator points of crop diversity change through time for each country (displayed visually in Fig. 1). Indicator 1 reflects the onset of diversification in each country, and was calculated as Breakpoint 1 (ψ1) in Eq. ( 2); this indicator therefore corresponds to the year in which notable changes in species richness began. Indicator 2 reflects the duration of the crop diversification period in each country, and was calculated as the difference between breakpoints 2 and 1 (i.e., ψ2-ψ1 from Eq. 2); this indicator therefore represents the duration of the period when crop prominent changes in crop diversity occurred. Finally, Indicator 3 reflects the rate at which crop diversity changed throughout the diversification period in each country; this indicator was calculated as the rate of crop diversity change (between ψ1 and ψ2), which in our models corresponded to the sum of the slopes (1) prior to the first breakpoint, and (2) between the first and second breakpoints (i.e., corresponding to b + c in Eq. 2). For each indicator we then calculated summary statistics as either mean ± standard deviations or median ± median absolute deviations (m.a.d.), where data was normally or log-normally distributed, respectively. Country values for each indicator were mapped using the mapCountryData function in the rworldmap R package [START_REF] South | rworldmap: A new R PACKAGE FOR MAPPING GLOBAL DATa[END_REF] .

Changes in crop evenness over time. Evaluations of temporal changes in crop evenness at national scales followed this same analytical approach as above. First, for each country-by-year combination we calculated Pielou's evenness index (J′)-which ranges from 0 to 1, with values closer to 0 indicating less evenness or greater abundance of a few dominant crop groups, and values closer to 1 representing more equitable abundances of crop groups-as:

where S is again crop richness, and H′ is the Shannon-Weiner diversity index calculated as:

(1) S = a + b × year (2) S = a + b year + c(year -ψ 1 × I year > ψ 1 + d year -ψ 2 × I year > ψ 2 (3) J ′ = H ′ ln (S)
where p i represents the relative proportion of the ith crop group for a given country-by-year combination. In these evenness calculations, all values of p i were estimated as the relative proportion of agricultural area (measured in ha) occupied by a given crop commodity group, within a country at a given year; this analytical approach was employed by Martin et al. [START_REF] Martin | Regional and global shifts in crop diversity through the Anthropocene[END_REF] when assessing crop group composition at supra-national scales. We then evaluated how J′ values changed in each country through time by replicating our stepwise modelling analyses above, substituting J′ for S in Eqs. ( 1) and ( 2), and extracting the same model indicators (Fig. 1). Finally, we calculated summary statistics and mapped each of these indicators, as described above.

Changes in crop composition across countries and over time. We used multivariate analyses to evaluate how temporal changes in S and J′ influenced crop composition across countries and over time. To do so, we created a community composition matrix whereby national-level crop assemblages were estimated for each of the country-by-year combinations. In this matrix, area harvested was taken as an approximation of the abundance of each crop group within each country-by-year combination (again following Martin et al. [START_REF] Martin | Regional and global shifts in crop diversity through the Anthropocene[END_REF] ). Since these abundances (or area harvested) across country-by-year combinations varied over orders of magnitude, we used non-metric multidimensional scaling (NMDS) to analyze and visualize spatial (country) and temporal (year) differences in crop diversity. Specifically, we used the vegan R package [START_REF] Oksanen | Community Ecology Package. R package version 2[END_REF] to calculate all 58,899,231 Bray-Curtis dissimilarities among all 10,854 data points (i.e., crop group composition in every country-by-year data point), as:

where BC jk represents the dissimilarity between the jth and kth community, x ij represents the abundance (i.e., area harvested) of crop group i in sample j, and x ik represents the abundance of crop group i in sample k. We then used a multivariate analysis of variance (i.e., an Adonis test), to test for significant differences in Bray-Curtis distances as a function of country, year, and a country-by-year interaction. Significance was assessed using a permutation test, with 99 permutations used.

Latitudinal gradients in crop richness. To test our hypotheses surrounding the presence of, and temporal changes in, latitudinal gradients in crop group diversity, we focused on 164 countries for which crop group diversity was available in both 1961 and 2017. For each of these two datasets, we fit a separate linear regression model that predicts crop group richness as a function of latitude (expressed as an absolute value) and a 2ndorder polynomial term for the 'latitude 2 ' variable. From both of these models, we extracted and compared latitude value at which crop group richness was estimated/ modelled to peak.

Predictors of change in crop diversity and composition.

We tested if Human Development Index (HDI) was correlated with patterns of change in crop diversity and composition. Briefly, the HDI is a composite index of four metrics related to socio-economic status, including life expectancy at birth, expected years of schooling for children at a school-centring age, mean years of schooling for adults ≥ 25 years of age, and logtransformed gross national income per capita. These values are then aggregated on a per country basis, into an HDI index that ranges from 0-1 with higher scores denoting higher performance in these indicators. We employed 2017 HDI values in our analysis here, in order to include the most countries possible in each analysis (since earlier HDI scores are less readily available) [START_REF]nlme: Linear and Nonlinear Mixed Effects Models. R package version[END_REF] .

We then used linear mixed effects models to test if patterns of change in crop diversity and evenness varied systematically with HDI values. This entailed fitting six linear mixed models, where each of our six indicators (i.e., Indicators 1-3 for both S and J′) were predicted as a function of HDI; these models also accounted for potential spatial autocorrelation in Indicator values by including the FAO-defined continent identity and FAOdefined region identity of each country, as a nested random variable. Models were fit using the lme function in the nlme R package [START_REF]nlme: Linear and Nonlinear Mixed Effects Models. R package version[END_REF] . We then estimated the proportion of variation in each indicator that is explained by HDI, continent identity, and region identity, using the varcomp function in the ape R package [START_REF] Paradis | Analyses of phylogenetics and evolution in R language[END_REF] -which partitioned explained variation across continents and regions-as well as the sem.model.fits function in the piecewiseSEM R package [START_REF] Lefcheck | piecewiseSEM: Piecewise structural equation modelling in R for ecology, evolution, and systematics[END_REF] -which partitioned explained variation across the fixed (i.e., model intercept and HDI) vs. random (i.e., continent and region) effects. Due to differences in HDI data availability and in the number of piecewise models that converged, n = 152 countries for all models of S indicators and n = 139 countries for all models of J′ indicators. Log-transformed values of Indicators were used in these analyses where they better approximated a log-normal distribution, as determined using the fitdistrplus function in the fitdistrplus R package [START_REF] Delignette-Muller | fitdistrplus: An R package for fitting distributions[END_REF] .

Figure 1 .

 1 Figure 1. Schematic representation of three indicators of change in crop commodity group diversity, derived from piecewise models predicting crop commodity group richness as a function of year. Detailed explanations of Indicator 1-3 are presented in the Methods section. Data shown here as the example is from Canada, with black dots representing the number of commodities reported by the Food and Agricultural Organization, for a given year. Black trendline represents the piecewise model fit, gray bands represent the 95% confidence limits surrounding the model, and red lines represent model parameters and indicators derived from the model. Note: the figure presented here demonstrates changes in crop commodity group richness (S), though this framework was also employed for assessing change in crop group evenness (J′).

Figure 2 .

 2 Figure 2. Maps and histograms of three indicators of crop commodity group richness (S) change across 165 countries.Values for all three indicators for each country were derived from piecewise linear models predicting S as a function of year (see Fig.1for example). Countries coloured gray in the maps were those where either data was not available or the piecewise models failed to converge (denoted in TableS1). Histograms and associated descriptive statistics for each indicator are also presented, with means (± s.d.) or medians (± m.a.d.) denoted visually by the points and error bars below the histograms. All piecewise model parameters for each country are presented in TableS1. Maps were generated using the mapCountryData function in the rworldmap R package[START_REF] South | rworldmap: A new R package for mapping global data[END_REF] .

Figure 4 .

 4 Figure 4. Maps and histograms of three indicators of crop commodity group evenness (Pielou's evenness index (J′)) across 185 countries.Values for all three indicators for each country were derived from piecewise linear models predicting J′ as a function of year, where harvested area (in ha) was used to approximate group abundance. Countries coloured gray in the maps were those where either data was not available or the piecewise models failed to converge (denoted in TableS2). Histograms and associated descriptive statistics for each indicator are also presented, with means (± s.d.) or medians (± m.a.d.) denoted visually by the points and error bars below the histograms. All piecewise model parameters for each country are presented in TableS2. Maps were generated using the mapCountryData function in the rworldmap R package[START_REF] South | rworldmap: A new R package for mapping global data[END_REF] .

Intercept Latitude Latitude 2 Model r 2 (p value) Latitude with modelled peak richness

	1961	16.1 (3.9) 1.3 (0.3)	-0.02 (0.006) 0.105 (0.0001)	37.04°2
	017	23.6 (4.5) 1.7 (0.4)	-0.02 (0.006) 0.129 (< 0.0001)	36.39°F

igure 3. Latitudinal patterns in crop group richness across 164 countries in 1961 and 2017.

Table 2 .

 2 Results from an Adonis test evaluating changes in crop commodity group composition at national scales from 1961 to 2017. The distance matrix employed in this analysis was based on non-metric multidimensional scaling, whereby cultivation area was used as an estimate for crop group abundance. Results are presented visually in FigureS1.

	Parameter	D.F	Sum of squares	Mean squares	Model F	r 2	p value
	Year	1	7.2	7.3	410.8	0.002	0.01
	Country	175	3862.8	22.1	1250.8	0.895	0.01
	Year*country	175	273.2	1.6	88.5	0.063	0.01
	Residuals	9680	170.8	0.02	-	0.04	-
	Total	10,031	4314.2	-	-	1	-

Table 3 .

 3 Results

	Variables		Fixed factors		Variance components (proportion explained)
	Metric	Indicator Intercept	HDI	Fixed effects Continent Region Unexplained
		1	1982.8 (5.0) 2.0 (6.5)	0.001	0.101	0.074	0.825
	Crop group richness (S)	log-2	1.8 (0.5)	0.4 (0.6)	0.003	0.015	0.000	0.982
		log-3 + 10 2.8 (0.1)	-0.5 (0.2)	0.068	0.074	0.000	0.925
		log-1	7.6 (0.003)	-0.004 (0.004) 0.009	0.234	0.04	0.716
	Crop group evenness (J′)	log-2	2.4 (0.5)	0.2 (0.6)	0.001	0.138	0.007	0.855
		3	0.001 (0.01)	-0.003 (0.01)	0.001	0.014	0.000	0.986

and variance components predicting Indicators of crop group diversity (S) and evenness (J′) change, as a function of Human Development Index (HDI), continent identity, and region identity. In these models an intercept and slope related to HDI were included as fixed effects, while region within continent were included as nested random effects. Statistically significant (p < 0.05) fixed effect model parameters are highlighted in bold, and data transformations (based on distribution fitting) precede the Indicator numbers (as per Figs.

2 and 4

).

Scientific Reports | (2021) 11:20361 | https://doi.org/10.1038/s41598-021-99728-5

Scientific Reports | (2021) 11:20361 | https://doi.org/10.1038/s41598-021-99728-5

Acknowledgements

This research was undertaken in part based on funding from Natural Sciences and Engineering Research Council of Canada (NSERC) Discovery Grants to Adam Martin, Marney Isaac, and Marc Cadotte. MEI is also supported by supported by NSERC's Canada Research Chairs program, while MWC is supported by the TD Professor of Urban Forest Conservation and Biology Chair program, the Canada Foundation for Innovation, and the Ontario Research Fund. Denis Vile is supported by the European Union and the Region Languedoc-Roussillon "Chercheur d' Avenir" (FEDER FSE IEJ 2014-2020; Grant Project "APSEVIR"). Cyrille Violle was supported by the European Research Council (ERC) Starting Grant Project StG-2014-639706-CONSTRAINTS. Funding from the University of Toronto Scarborough's (UTSC) Vice Principal Research Office in support of UTSC's Centre for Agroecosystems Research also supported this research.

Author contributions

R.O.M. lead data consolidation, analysis, and visualization, and wrote the manuscript. A.R.M. conceived the study, contributed to data analysis and visualization, and co-wrote the manuscript. M.W.C, M.E.I., R.M., D.V., and C.V. contributed to manuscript editing and data analysis.

Competing interests

The authors declare no competing interests.


Yuta Shimizu

Dr Koichiro Suekuni
email: suekuni.koichiro.063@m.kyushu-u.ac.jp

Dr Hikaru Saito

Dr Pierric Lemoine

Dr Emmanuel Guilmeau

Bernard Raveau

Dr Raju Chetty

Dr Michihiro Ohta

Toshiro Takabatake

Michitaka Ohtaki

Prof M Ohtaki

Prof. B Raveau

Prof T Takabatake

Synergistic Effect of Chemical Substitution and Insertion on the Thermoelectric Performance of Cu26V2Ge6S32 Colusite

Keywords: thermoelectric, sulfide, colusite, elemental substitution, interstitials, thermal conductivity

Copper-based sulfides are promising materials for thermoelectric applications, which can convert waste heat into electricity. This study reports the enhanced thermoelectric performance of Cu26V2Ge6S32 colusite via substitution of antimony (Sb) for germanium (Ge) and introduction of copper (Cu) as an interstitial atom. The crystal structure of the solid solutions and Cu-rich compounds were analyzed by powder X-ray diffraction and scanning transmission electron microscopy. Both chemical approaches decrease the hole carrier concentration, which leads to a reduction in the electronic thermal conductivity while keeping the thermoelectric power factor at a high value. Furthermore, the interstitial Cu atoms act as phonon scatterers, thereby decreasing the lattice thermal conductivity. The combined effects increase the dimensionless thermoelectric figure of merit ZT from 0.3 (Cu26V2Ge6S32) to 0.8 (Cu29V2Ge5SbS32) at 673 K.

Introduction

Renewable energy sources, which include sunlight, heat, and vibration, play a crucial role in realizing a sustainable society. One of the promising applications for harvesting heat to generate electricity is the thermoelectric (TE) generator, which is composed of solid-state elements. [START_REF] Bell | [END_REF][2][3] The conversion from thermal to electrical energy in solids is based on the Seebeck effect. The electronic potential difference (electromotive force) ΔV is generated in proportion to the temperature difference ΔT (ΔV = SΔT, where S is the Seebeck coefficient) between two ends of an element. The electromotive force is utilized to supply electricity to an electronic device. Therefore, highly efficient TE conversion requires an element material to have large S. Other requisite characteristics of a TE material include low electrical resistivity ρ to reduce the internal resistance of the generator and low thermal conductivity κ to maintain sufficient ΔT. By combining these parameters (ρ, S, κ), the performance of a TE material can be determined using the expression: ZT = S 2 Tρ -1 κ -1 , where ZT is referred to as the dimensionless figure of merit, T is the absolute temperature, and κ is the sum of its electronic component κele and lattice component κlat. In addition to high performance (high ZT), naturalabundance and non-toxic characteristics are desired for TE materials used in large-scale applications. As such TE materials, Cu-S-based (degenerate) semiconductors have emerged in the recent decade. [4,5] Examples include Cu2ZnSnS4, [6] CuFeS2, [7,8] Cu1.8S, [9] Cu2S, [10] Cu12Sb4S13, [11][12][13] Cu3SbS4, [14] Cu5FeS4, [15][16][17][18] Cu2SnS3, [19] Cu4Sn7S16, [20] Cu8Fe3Sn2S12, [21] CuFe2S3, [22] Cu22Fe8Ge4S32, [23] Cu3PS4, [24] and Cu26V2Ge6S32. [25,26] The last example is a synthetic member of the colusite family: Cu26T2M6S32 (T = Ti, V, Nb, Ta, Cr, Mo, W; M = Ge, Sn, Sb). [25][26][27][28][29][30] Colusites have attracted increasing attention due to their high ZT value ranging from 0.5-1.0 at 673 K. [26][27][28][29][30][31][32][33][34][35][36][37][38][39][40][41][42][43] Cu26T2M6S32 colusites crystallize in a cubic structure (P4 � 3n), which is composed of three types of tetrahedral units: CuS4, TS4, and MS4. [44,45] The CuS4 and MS4 tetrahedra share corners to construct a three-dimensional framework structure, which is completed by the TS4 tetrahedra sharing edges with CuS4 (Figure 1). The formal valence of the atoms can be represented as Cu 32 for T = Cr, Mo, and W. The non-magnetic and p-type degenerate semiconducting characteristics of the colusites with T = V, Cr, Mo, and W were confirmed by TE and magnetic susceptibility measurements, [25] Mössbauer spectroscopy, [46] and firstprinciples electronic-structure calculations. [26,28,36] These characteristics indicate the itinerant nature of 3d holes of Cu 2+ . [28,47] The calculations showed that the Fermi level (EF) lies at the top of the valence band (VB), which primarily comprises the hybridized orbitals of Cu-3d and S-3p. [26,28,36] The nature of the VB indicates that the Cu-S-based tetrahedral framework governs the electronic properties. Depending on the compositions, a high power factor S 2 ρ -1 of 0.5-1.9 mW K -2 m -1 at 673 K can be achieved in colusites. [26][27][28][29][30][31][32][33][34][35][36][37][38][39][40][41][42][43] The tetrahedral framework also causes low κlat ≤ 1 W K -1 m -1 due to its complexity and low-energy (~10 meV) optical phonon modes involving motions of Cu and S. [28,29,36] Further decrease in κlat can be achieved by introducing various types of defects into the crystal structure: interstitial defects, anti-site defects between cations, and splitted Cu site. [36,37,46] These defects are simultaneously produced by an excessive amount of sulfur sublimation upon heating above 973 K. [36,37] A good approach to introduce phonon scatterers (defects) into the colusite structure is by increasing the sintering temperature; [36,37] however, the difficulty in controlling sulfur loss prevents fine-tuning of the hole carrier concentration n. Therefore, an alternative strategy to reduce κlat is required. Generally, substitutional doping of a heavier atom for a constituent atom results in a decrease in κlat. [48] Another possible way to reduce κlat in colusites is to introduce interstitial defects (atoms), according to the aforementioned results. [37] However, the effect of interstitial atoms on κlat has not been sorely examined for colusites to date. Following these circumstances, this study aimed to clarify the effects of substituting Sb for Ge and Cu insertion on κlat in Cu26V2Ge6S32 colusite. Hence, novel Sb-substituted (Cu26V2Ge6-xSbxS32),

Cu-rich (Cu26+yV2Ge6S32), and combined Sb-substituted and Cu-rich (Cu26+yV2Ge6-xSbxS32 with x = 1) series were synthesized and their TE properties were investigated. Notably, these samples were sintered at 873 K, which is mandatory for limiting the sublimation of sulfur. [36,37,49] S3.

Results and Discussion

The powder X-ray diffraction (PXRD) patterns of the samples in different series of compositions, i.e., Cu26V2Ge6-xSbxS32 (x = 0-4), Cu26+yV2Ge6S32 (y = 0-4), and Cu26+yV2Ge6-xSbxS32 (x = 1, y = 0-3) are presented in Figure 2a and Figures S1-S3. For all samples, the positions and relative intensities of the diffraction peaks reasonably agree with that of a simulated pattern based on the colusite structure (Cu26V2Ge6S32, 𝑃𝑃4 � 3𝑛𝑛). Note that the systematic presence of weak-intensity diffraction peaks (at 26.4° and 43.7° for x = 0) is originated from CuKβ radiation passing through a Ni filter. The diffraction patterns for x = y = 0 (pristine sample) and y ≥ 1 (Cu-rich samples) exhibited small peaks from a secondary phase Cu2S, and that for x = 4 (Sb-substituted sample) showed a small peak from an unidentified phase (Figures S1-S3). For the Cu-rich series, the peak width became broader as y increased (Figures S2 andS3).

Rietveld refinement of the PXRD patterns of the Sb-substituted series (Cu26V2Ge6-xSbxS32) (Figure S4-S8 and Table S1), leading to low-reliability factors, confirmed the validity of the structural model of colusites, which includes seven crystallographic sites: three for Cu (6d, 8e, 12f), one for V (2a), one for Ge/Sb (6c) and two for S (8e, 24i). The lattice parameter a linearly increased from 10.5684(2) Å for x = 0 by 1.2% as x increased to 4, as shown in Figure 2b and Table S1. The lattice expansion can be attributed to the larger ionic radius of Sb 5+ (0.565 Å) [50] than that for Ge 4+ (0.390 Å). [51] In fact, the occupation factor of Sb at Ge (6c) site obtained from Rietveld refinement increased with x and reached 62% (i.e. Cu26V2Ge2.3Sb3.7S32) for x = 4 (Table S1). The increase in the Sb content was also clearly observed from the evolution of h + k + l = 2n + 1 peaks (e.g. 210, 320, see Figure 2a and Figure S1) [29] . It was further confirmed that this substitution increases the M(6c)-S(24i) (M = Ge, Sb) distance from 2.254 Å for x = 0 to 2.383 Å for x = 4 while keeping the average Cu-S distance intact (2.291-2.300 Å) (Table S1). This indicates the limited modification of Cu-S-based tetrahedral framework. It should be noted that the presence of Sb 5+ in tetrahedral coordination of sulfides was only reported for famatinite Cu3SbS4 [52] and the colusite sister compound Cu26Ti2Sb6S32. [29] For the Cu-rich series (Cu26+yV2Ge6S32 and Cu26+yV2Ge6-xSbxS32 with x = 1),

Rietveld refinement of the PXRD patterns confirmed a lattice expansion attributed to Cu insertion (Figure 2b). More specifically, a increased slightly for y = 1 and increased by 0.44% and 0.47% for y = 3 in Cu26+yV2Ge6S32 and Cu26+yV2Ge6-xSbxS32 (x = 1), respectively.

Saturation of the a value for Cu30V2Ge6S32 (y = 4) suggests a lower (actual) Cu content, which is consistent with the increase in the amount of the Cu-rich secondary phase (Cu2S) for y = 4

and reflects a limit of Cu insertion into the VGe-colusite structure. A possible mechanism of the lattice expansion for the Cu-rich series deals with the introduction of Cu atoms into interstitial sites, as previously reported for Cu26Nb2Sn6S32-δ colusite. [37] The fact that Cu insertion leads to electron doping (decrease in n), as shown below, strongly supports this mechanism. images show atomic columns. Along this direction, there are four types of columns labelled "4Cu", "1V", "3Cu+M" / "3Cu+Ge", and "4S" (see text). Here, the numeral denotes the number of atoms periodically contained in the columns, and M means "(4Ge + 2Sb) / 6".

Logarithmic intensities (levels of brightness) of columns along "lines" denoted by 1-5 in the expanded views are plotted in the graphs. The dashed lines in the graphs are eye guides for making peaks/humps clear. For (b) Cu26+yV2Ge6S32 (y = 3), the intensity along a line composed of the "3Cu + 1Ge" and "4S" columns shows a hump at an interstitial position.

Such a hump was hardly noticed for (a) Cu26V2Ge6S32 and (c) Cu26V2Ge6-xSbxS32 (x = 2).

To confirm the existence of interstitial Cu ions, scanning transmission electron microscopy (STEM) analysis was performed for pristine, Sb-substituted, and Cu-rich samples, as shown in Figure 3 [110] direction (Figure 3), there were four types of columns: pure copper columns labelled "4Cu", pure vanadium columns labelled "1V", pure sulfur columns labelled "4S", and mixed columns containing besides copper either mixtures of Ge and Sb labelled "3Cu + 1M" or only Ge labelled "3Cu + 1Ge". Here, the numeral denotes the number of atoms periodically contained in the columns and M means "(4Ge + 2Sb) / 6". Bearing in mind that the intensity of the spots corresponding to the columns is approximately proportional to the square of the atomic number Z (29 for Cu, 32 for Ge, 51 for Sb, 23 for V and 16 for S), one can expect a decrease in their relative intensity from "3Cu + 1M" to "3Cu + 1Ge" to "4Cu" to "4S" to "1V" columns and that the intensities of the "3Cu + 1Ge" and "4Cu" columns are comparable. As a result, the "3Cu + 1M" / "3Cu + 1Ge", "4Cu" and "4S" columns are clearly identified in the STEM images. As shown in the line profiles of the intensity (Figure 3 and Figures S10 andS12), the relative intensities of the "3Cu + 1M" column were higher than that of the "3Cu + 1Ge" column, which is consistent with the substitution of Sb for Ge determined from PXRD analyses. The line profile also allows us to see the "1V" columns, whose position agrees with that of the colusite structure.

Remarkably, for the Cu-rich sample (Cu26+yV2Ge6S32 with y = 3), the logarithmic intensity along the line composed of "3Cu + 1Ge" and "4S" columns showed a hump at an interstitial position for the colusite structure (Figure 3c). Such a hump was hardly noticed in the pristine and Sb-substituted samples, Cu26V2Ge6S32 and Cu26V2Ge6-xSbxS32 with x = 2 (Figures 3a,3b). These results strongly indicate the existence of atoms at interstitial sites in the Cu-rich sample. This position may be consistent with the presence of Cu at an interstitial 24i site (0.255, 0.215, 0.032), as determined by the powder synchrotron X-ray diffraction for Cu26Nb2Sn6S32-δ. [37] However, signals from interstitial atoms were observed at a part of "columns" of the interstitial sites. This suggests a non-uniform arrangement (or aggregation) of interstitial atoms, which might result in the broadening of diffraction peaks (Figure 2a and The existence of interstitial atoms was further confirmed by Rietveld refinement of the PXRD pattern for Cu26+yV2Ge6S32 with y = 3 (Figure S13). Refinement based on the colusite structure showed a residual electron density at an interstitial 24i site (0.236, 0.235, 0.007). Then, refinement based on a colusite structure with interstitial Cu atom at the 24i site (atomic coordinates being fixed) was performed, which significantly improved the reliability factors (Table S2). The refined occupation factor at the 24i site was 0.042 (1), which corresponds to an extra 1.0 Cu atom per formula unit, i.e. Cu27V2Ge6S32. Although the quality of the PXRD pattern was not sufficient for quantitative analysis (e.g. broadened peaks), the refinement confirmed partial occupation of a 24i site (~0.24, ~0.24, ~0.01) by extra atoms. Next, the TE properties for the three series of samples were discussed (except for x = 4 and y = 4). All samples had relative densities higher than 95%, as confirmed by scanning electron microscopy (SEM) (Figure S14). For the Sb-substituted series (Cu26V2Ge6-xSbxS32), all samples exhibited the p-type metallic behavior of degenerated semiconductors. The pristine sample (x = 0) had low S (18 μV K -1 ) and ρ (1.7 μΩ m) values at 300 K (Figures 4a and4b), which is in agreement with previous studies. [25,28] The values of S and ρ increase with

x and reaches 116 μV K -1 and 17 μΩ m, respectively, for x = 3. The increase in S and ρ results from a decrease in n from 6.8 × 10 21 cm -3 for x = 0 to 1.9 × 10 21 cm -3 for x = 3 (i.e. electron doping) (Table S3), which can be attributed to the larger valence of Sb 5+ compared to Ge 4+ .

As a result, the power factor S 2 ρ -1 increased from 1.2 mW K -2 m -1 for x = 0 to 1.4 mW K -2 m -1 for x = 1 and x = 2, and then decreased to 1.0 mW K -2 m -1 for x = 3 at 673 K (Figure 4c).

The increase in ρ led to a strong reduction in κele, that manifested itself in the form of a drastic decrease in κ (Figure 4d). Consequently, the combination of enhanced S 2 ρ -1 and reduced κ resulted in the enhancement of ZT from 0.3 (x = 0) to 0.7 (x = 3) at 673 K for the Sbsubstituted samples (Figure 4e). For the Cu-rich series (Cu26+yV2Ge6S32), the values of S and ρ increased with y (Figures 5a and5b), and it is consistent with a decrease in n from 6.8 × 10 21 cm -3 for y = 0 to 4.1 × 10 21 cm -3 for y = 3 at 300 K (Table S3). The decrease in n originates from that the insertion of Cu + doped electron into the VB. As a result, S 2 ρ -1 remains at a high level (1.2-1.4 mW K -2 m -1 ) at 673 K (Figure 5c). The value of κ decreases with y due to the reduction in κele (Figure 5d). Consequently, ZT values increase from 0.3 (y = 0) to 0.6 (y = 3) at 673 K (Figure 5e). For the combined Sb-substituted and Cu-rich series (Cu26+yV2Ge6-xSbxS32 with x = 1), the increase in S and ρ (Figure 6a and 6b) with increasing y is consistent with the aforementioned results for Cu26+yV2Ge6S32 (Figure 5a and5b). Remarkably, an increase in ρ was observed at temperatures below 400 K for Cu26+yV2Ge6-xSbxS32 (x = 1, y =3) (Figure 6b), unlike the metallic behaviour observed in ρ for Cu26V2Ge6-xSbxS32 (x = 3) (Figure 4b).

Because the samples of x = 1, y = 3 and x = 3 have similar n at 300 K (i.e., 1.7 × 10 21 cm -3 for

x =1, y = 3; 1.9 × 10 21 cm -3 for x = 3), the enlargement of ρ for x =1, y = 3 can be attributed to a decrease in the electron mobility µH from 1.9 cm 2 V -1 s -1 to 1.0 cm 2 V -1 s -1 (Table S3).

Also, because the electronic properties of colusite are governed by the Cu-S tetrahedral framework, [39] the presence of interstitial atoms probably affects µH in the conductive network.

A similar upturn in ρ was observed in cation-rich [Cu26Cr2Ge6]1+δS32 colusites. [40] As a result of the aforementioned variation of y for Cu26+yV2Ge6-xSbxS32 (x = 1), the value of S 2 ρ -1 at 300 K remains unchanged, whereas it increases at 673 K for y = 1 and decreases for y ≥ 2 (Figure 6c). The sample produced at x = 1, y = 1 showed the highest S 2 ρ -1 value of 1.5 mW K -2 m -1 at 673 K compared to the other samples in this study. The value of κ (and κele) decreased as y increased (Figure 6d). As a result, the reduction in κ boosts ZT at 673 K from 0.4 (x = 1, y = 0) to 0.8 (x = 1, y = 3) (Figure 6e). The ZT value of the samples produced at x = 1, y = 3 in the Cu26+yV2Ge6-xSbxS32 series was higher than that of x = 3 in the Cu26V2Ge6-xSbxS32 series (Figure 4e) despite similar values of S and ρ at 673 K (Figures 4 and6). This comparison indicates that the reduction in κlat via Cu insertion contributes to the enhancement of ZT. To determine the effects of Sb substitution and Cu insertion on κlat, the value of κlat was estimated by subtracting κele = LTρ -1 from κ, where L = 1.5 + exp (-|S|/116). [53] The equation regarding L is based on the acoustic phonon scattering of electrons, which is consistent with the metallic increase in ρ(T) at temperatures between 400 K and 673 K for all the samples (Figures 4b, 5b κlat slightly increased for y ≥ 1, whereas for the Cu26+yV2Ge6-xSbxS32 (x = 1) series, κlat decreased for y ≥ 1. These opposite trends for the Cu-rich series imply that the variation of κlat is not influenced by the existence of Cu2S (secondary phase), whose amount increased with y in both series. The variations of κlat can originate from different effects/mechanisms: (i) the substitutional doping of Sb for Ge, (ii) the interstitial doping of Cu, (iii) uncertainty in the estimation of L, due to the complexity of the actual electronic structure [28] and (iv) phonon scattering by electrons. [54] Hence, we chose three pairs of samples in the Cu26+yV2Ge6-xSbxS32 series with similar S and ρ values (equivalent n and EF) to minimise the influence of (iii) and

(iv) on κ. Each pair includes Sb-substituted/non-Cu-rich and Cu-rich samples (i.e. x = 1, y = 0 (n = 5.5 × 10 21 cm -3 at 300 K) vs. x = 0, y = 2 (5.1 × 10 21 cm -3 ); x = 2, y = 0 (4.0 × 10 21 cm -3 ) vs. x = 0, y = 3 (4.1 × 10 21 cm -3 ); x = 3, y = 0 (1.9 × 10 21 cm -3 ) vs. x = 1, y = 3 (1.7 × 10 21 cm -3 )). As shown in Figure 7, although the substitution of Sb for Ge has a negligible effect on κlat, all Cu-rich samples showed lower κlat compared with the non-Cu-rich samples. This result demonstrates that the interstitial Cu ions act as phonon scatterers in the crystal structure, leading to the reduction of κlat. Furthermore, the decrease in κlat at approximately 300 K was more pronounced in the samples showing larger upturn in ρ(T). This means that the insertion of Cu has a significant impact on both enhanced phonon scattering and reduced µH (Table S3).

Further investigations (e.g., Mӧssbauer spectroscopy, extended X-ray absorption fine structure spectroscopy) are required to investigate the local structure modification of the colusite structure.

Conclusion

In this work, we performed the substitution of Sb for Ge and the insertion of Cu into the crystal structure of Cu26+yV2Ge6-xSbxS32 colusite. Both chemical approaches led to a decrease in n and κele, while maintaining S 2 ρ -1 ≥ 0.9 mW K -2 m -1 (at x ≤ 3, y ≤ 3).

Furthermore, the interstitial Cu atoms act as phonon scatterers, leading to the suppression of κlat. These combined effects boost ZT, which reaches a value of 0.8 at 673 K for x = 1, y = 3 sample. Thus, the Cu insertion into corner sharing three-dimensional tetrahedral framework structures of Cu-S-based materials is a powerful chemical approach for improving TE performance via tuning of n and reducing of κlat.

Experimental Section

The samples of Cu26V2Ge6-xSbxS32 (x = 0-3), Cu26+yV2Ge6S32 (y = 0-4) and Cu26+yV2Ge6-xSbxS32 (x = 1, y = 0-3) were prepared by directly reacting the constituent elements at 1373 K followed by heat treatment at 873 K as reported for the samples of Cu26Nb2Sn6S32. [37] The obtained samples were pulverized using a planetary ball mill PXRD data were collected in the range of 10° ≤ 2θ ≤ 100° using an X-ray diffractometer (MiniFlex600, RIGAKU) with a Cu Kα radiation source. PXRD data analyzes were performed by Rietveld refinement using the FullProf and WinPlotr software packages. [55,56] Zero-point shift, lattice parameter, peak shape parameters, asymmetry parameters, fractional atomic coordinates and isotropic displacement parameters were refined after manually removing the background contribution. The density of the sample was checked by SEM using a microscope (JCM-6000Plus NeoScope, JEOL). Atomic-resolution observations were performed using a transmission electron microscope (Titan Cubed 60-300 G2, Thermo

Fisher Scientific) which is equipped with a spherical aberration corrector (DCOR, CEOS) for the probe-forming lens system. The microscope was operated in the STEM mode at an accelerating voltage of 300 kV. The convergence semi-angle of the electron probe was set to 18 mrad. The typical probe diameter was less than 0.1 nm. An ADF detector was positioned to detect scattered electrons with an angular range from 38 to 184 mrad. Each sintered bulk sample was crushed and dispersed to ethanol in Ar atmosphere to avoid oxidation. Then, it was dropped on a carbon supporting film and was immediately inserted in the electron microscope.

S and ρ were simultaneously measured at T between 300 K and 673 K in a measurement system (ZEM-3, ADVANCE RIKO) in a helium atmosphere. Thermal diffusivity α and specific heat CP at 300-673 K were measured simultaneously in a measurement system (LFA 457 MicroFlash, Netzsch) using the laser-flash method in a flowing Ar atmosphere. In this study, the sample was coated with a thin layer of graphite. The absolute values of CP were derived from the comparison between the measured values and CP values measured for a standard sample of Pyroceram 9606 (Netzsch) (Figure S16). These data were used to calculate κ as αCPds, where ds represented sample density estimated from the dimensions and weight of the sample. Note that the obtained Cp values reasonably agreed among samples, which is comparable to the Dulong-Petit value at 300 K and higher by ≤ 16% at 673 K (Figure S16). The Hall-effect measurements were performed using a four-probe DC method on a laboratory-built system, with a permanent magnet generating a magnetic field of Table S1. Results from the Rietveld refinements of power X-ray diffraction patterns for Cu26V2Ge6-xSbxS32 samples (space group 𝑃𝑃4 � 3𝑛𝑛). Table S2. Results from Rietveld refinements of powder X-ray diffraction patterns for y = 0 and y = 3 samples of Cu26+yV2Ge6S32 (space group 𝑃𝑃4 � 3𝑛𝑛).

Cu26V2Ge6-xSbxS32 x = 0 x = 1 x = 2 x = 3 x = 4 a (Å) 10
Cu26+yV2Ge6S32 y = 0 y = 3 y = 3 a (Å)

Figure 1 .

 1 Figure 1. Crystal structure of Cu26T2M6S32 (T = Ti, V, Nb, Ta, Cr, Mo, W; M = Ge, Sn, Sb)

Figure 2 .

 2 Figure 2. (a) Powder X-ray diffraction patterns and (b) lattice parameters a for

Figure 3 .

 3 Figure 3. Annular dark-field scanning transmission electron microscopy (STEM) images for

  and Figures S9-S12. The annular dark-field (ADF) STEM images for Cu26V2Ge6S32, Cu26V2Ge6-xSbxS32 with x = 2, and Cu26+yV2Ge6S32 with y = 3 showed wellaligned atom columns with intensities (levels of brightness in the image) that are mostly consistent with those expected from the pristine colusite (Cu26V2Ge6S32) structure. Along the

Figure

  Figure S2). Signals from interstitial ions in the STEM images along the [100] and[111]

Figure 4 .

 4 Figure 4. (a) Seebeck coefficient S, (b) electrical resistivity ρ, (c) power factor S 2 ρ -1 , (d)

Figure 5 .

 5 Figure 5. (a) Seebeck coefficient S, (b) electrical resistivity ρ, (c) power factor S 2 ρ -1 , (d)

Figure 6 .

 6 Figure 6. (a) Seebeck coefficient S, (b) electrical resistivity ρ, (c) power factor S 2 ρ -1 , (d)

Figure 7 .

 7 Figure 7. (a) Seebeck coefficient S, (b) electrical resistivity ρ and (c) lattice thermal

  , and 6b). Conversely, single parabolic band assumed for the equation may be too simple to describe the VB of colusites. As shown in Figure S15, the obtained κlat for Cu26V2Ge6-xSbxS32, Cu26+yV2Ge6S32, and Cu26+yV2Ge6-xSbxS32 (x = 1) exhibited complicated variations with different x and y values. More specifically, κlat for Cu26V2Ge6-xSbxS32 increased for x = 2 and decreased for x = 3. For the Cu26+yV2Ge6S32 series,

(Pulverisette 7

 7 premium line, Fritsch) operating at room temperature at a disk rotation speed of 450 rpm for 1 h. The sample powder was put into a jar together with seven balls of 10 mm diameter in an Ar atmosphere. The jar and balls were made of tungsten carbide (WC). The pulverized sample was loaded into a WC die with an inner diameter of 10 mm. The sample was sintered at 823/873 K for 40 min in a flowing N2 atmosphere under a uniaxial pressure of 200 MPa in a hot-press sintering furnace (PLASMAN CSP-I-03121, S. S. Alloy). The relative density of the samples was higher than 95% of the theoretical density. The sintered sample was cut and polished into bars and discs to measure the TE properties.

  Figure S2. Powder X-ray diffraction patterns for Cu26+yV2Ge6S32 samples. Arrows denote peaks from Cu2S. The simulated pattern based on the colusite (Cu26V2Ge6S32) structure is shown at the bottom. The right panel shows expanded views of 622 peaks.

Figure S3 .

 S3 Figure S3. Powder X-ray diffraction patterns for Cu26+yV2Ge6-xSbxS32 (x = 1) samples. Arrows denote peaks from Cu2S. The simulated pattern based on the colusite (Cu26V2Ge6S32) structure is shown at the bottom. The right panel shows expanded views of 622 peaks.

Figure S4 .

 S4 Figure S4. Rietveld refinement of power X-ray diffraction patterns for Cu26V2Ge6S32.

Figure S5 .

 S5 Figure S5. Rietveld refinement of power X-ray diffraction patterns for Cu26V2Ge6-xSbxS32 with x = 1.

Figure S6 .

 S6 Figure S6. Rietveld refinement of power X-ray diffraction patterns for Cu26V2Ge6-xSbxS32 with x = 2.

Figure S7 .

 S7 Figure S7. Rietveld refinement of power X-ray diffraction patterns for Cu26V2Ge6-xSbxS32 with x = 3.

Figure S8 .

 S8 Figure S8. Rietveld refinement of power X-ray diffraction patterns for Cu26V2Ge6-xSbxS32 with x = 4.

Figure S9 .

 S9 Figure S9. Annular dark-field scanning transmission electron microscopy images along the (a) 100, (b) 110 and (c) 111 directions for Cu26V2Ge6-xSbxS32 (x = 2) sample.

Figure S10 .

 S10 Figure S10. Line profiles of intensity for atomic columns along the (a) 100, (b) 110 and (c) 111 directions depicted in Figure S9 and the corresponding view of the crystal structure.

Figure S11 .

 S11 Figure S11. Annular dark-field scanning transmission electron microscopy images along the (a) 100, (b) 110 and (c) 111 directions for Cu26+yV2Ge6S32 (y = 3) sample.

Figure S12 .

 S12 Figure S12. Line profiles of intensity for atomic columns along the (a) 100, (b) 110 and (c) 111 directions depicted in Figure S11 and the corresponding view of the crystal structure.

Figure S13 .

 S13 Figure S13. Rietveld refinement of power X-ray diffraction patterns for Cu26+yV2Ge6S32 with y = 3.

Figure S14 . 1 xFigure S15 .

 S141S15 Figure S14. Scanning electron microscopy images for Cu26V2Ge6-xSbxS32, Cu26+yV2Ge6S32, and Cu26+yV2Ge6-xSbxS32 (x = 1) samples.

Figure S16 .

 S16 Figure S16. Specific heat Cp in units of J g -1 K -1 and J mol -1 K -1 for Cu26V2Ge6-xSbxS32, Cu26+yV2Ge6S32, and Cu26+yV2Ge6-xSbxS32 (x = 1) samples.

  Rexp 4.09 ; 4.96 4.19 ; 5.23 3.98 ; 5.25 4.11 ; 5.49 4.62 ; 5.60 RBragg ; RF 2.62 ; 7.54 3.05 ; 6.45 2.20 ; 5.37 2.72 ; 5.17 3.03 ; 6.29

		.5684(2) 10.6026(2) 10.6334(2) 10.6671(2) 10.6988(2)
	SOF_Ge(6c)	1.00	0.87(2)	0.69(2)	0.50(2)	0.38(2)
	SOF_Sb(6c)	-	0.13(2)	0.31(2)	0.50(2)	0.62(2)
	x_Cu(8e)	0.248(1)	0.248(1)	0.247(1)	0.245(1)	0.246(1)
	x_Cu(12f)	0.260(1)	0.260(1)	0.257(1)	0.257(1)	0.256(1)
	x_S(8e)	0.126(1)	0.127(1)	0.125(1)	0.123(1)	0.124(1)
	x_S(24i)	0.371(1)	0.371(1)	0.372(1)	0.375(1)	0.376(1)
	y_S(24i)	0.373(1)	0.373(1)	0.372(1)	0.370(1)	0.368(1)
	z_S(24i)	0.121(1)	0.122(1)	0.122(1)	0.126(1)	0.128(1)
	χ 2	0.678	0.641	0.575	0.560	0.680
	Rwp ; Distances (Å)
	V(2a)-S(8e) × 4	2.308	2.325	2.297	2.267	2.296

Acknowledgements

The authors would like to thank Enago (www.enago.jp) for the English language review. This work was financially supported by JSPS KAKENHI Grant No. JP20H02440 (K. S.), and grants from the International Joint Research Program for Innovative Energy Technology funded by METI and Research and Development Program for Promoting Innovative Clean Energy Technologies Through International Collaboration funded by NEDO. E.G. and P.L. acknowledge the financial support of CNRS through the International Emerging Actions program (EXPRESS project).

Supporting Information

Supporting Information is available from the Wiley Online Library or from the author.