It is generally believed that marsupials are more primitive mammals than placentals and mainly solitary living, representing the ancestral form of social organization of all mammals. However, field studies have observed pair- and group-living in marsupial species, but no comparative study about their social evolution was ever done. Here we describe the results of primary literature research on marsupial social organization which indicate that most species can live in pairs or groups and many show intra-specific variation in social organization. Using Bayesian phylogenetic mixed-effects models with a moderate phylogenetic signal of 0.18 we found that solitary living is the most likely ancestral form (35% posterior probability), but has high uncertainty, and the combined probability of partly sociable marsupial ancestors (65%) should not be overlooked. For Australian marsupials, group-living species were less likely to be found in climates representing tropical rainforest, and species with a variable social organization were associated with low and unpredictable precipitation representing deserts. Our results suggest that modern marsupials are more sociable than previously believed and that there is no strong support that their ancestral state was strictly solitary living, such that the assumption of a solitary ancestral state of all mammals may also need reconsideration.

The ancestral state of the social organization of all marsupials: Our aim was to create a database based exclusively on published primary data, without consideration of interpretations or generalisations of the authors, and in doing so contribute to improving the quality of future comparative studies. We only considered primary literature that reported actual field data on social organization in their methods or results section. This was important for the main aim of our study which is to consider IVSO and not only the most common form of social organization believed to occur in one species. Reviews making generalisations were not considered.

We searched for publications about social organization of marsupials on Web of Science and Google Scholar from June 2020 until April 2021. Our search included all 345 marsupial species categorized by the IUCN (International Union for Conservation of Nature) database in 2021 (https://www.iucnredlist.org/). Each species was searched by its scientific name (genus and species) and the keyword “social” (e.g., Acrobates pygmaeus AND social). If no results were found, a second search was conducted using only the scientific name (genus and species). For each study, the title and abstract were read to determine whether the study was based on a wild population and if it might contain data about social organization. By reading the article titles from the search results we marked 697 articles, and downloaded 456 of them after reading the abstract. For 105 articles (mainly old articles from local journals) we were not able to obtain a PDF or copy; thus, we could not check them for suitability nor add them to our database.

To address our main interest in identifying primary data on IVSO (deviation from the main form of social organization), methods, results, figures and tables of all 456 articles were checked. Further, the full text was searched for the following keywords: "social", “solitary”, “pair”, “group”. 253 articles did not contain useable data on social organization. Based on our criteria, 83 articles could not be included in the analyses, as the authors only stated the main form of social organization, but did not present the data on the composition of social units, the sex of individuals, occupancy of sleeping sites, home range overlap or the proportion of the individuals marked and studied in the study area [1]. In other words, these studies were excluded because either IVSO was ignored as a possibility, or if the existence of IVSO was reported, it was impossible to determine the degree of it. Supplementary materials 8 presents an alternative analysis focusing only on the main form of social organization ignoring IVSO, including these studies excluded from our main analysis. As would be predicted, this model overemphasizes the probability of solitary living being ancestral.

We found data on social organization that matched our inclusion criteria in 120 of all articles. The data in these articles were recorded at the population level (N=149 populations) and covered 65 species. The phylogenetic distribution of those species shows that while there is very limited knowledge for the Ameridelphia, the available data for the Australidelphia are relatively evenly distributed among families.

We were able to classify the social organization of each social unit reported in the 120 articles as one of six possible forms: (1) solitary, (2) pair-living, and four forms of group-living: (3) single male multiple female group, (4) single female multiple male group, (5) sex-specific group (group of only males or only females), and (6) multi-male multi-female group. For analyses we used the category “stable group” if all social units showed the same category of group-living, while if more than one category occurred, this was categorised as IVSO (see details below). Social organization was only based on the number of adults present and we did not consider pups and juveniles (Joeys). Solitary individuals were recorded separately by sex; a social unit was recorded as solitary only when both sexes were solitary. Many species show sex-specific dispersal. Thus, when solitary living was only reported for individuals of one sex, this was not considered evidence for a solitary social organization, since the data might represent dispersing individuals. As most species have dispersing individuals that are solitary for a short time, including this transitional phase as a separate social organization would basically mean that all species show IVSO, in which case it would not make sense to study why and when it occurs. Instead, apart from studies on IVSO, separate studies on the proximate causes and ultimate function of dispersal are needed. To facilitate comparison to pairs (one male and one female), the number of solitary social units was determined by the sex with the smaller number of solitary individuals (e.g., when 10 solitary males and five solitary females were observed, we recorded five solitary social units). The same method was applied to sex-specific groups.

Intra-specific variation in social organization (IVSO) was recorded when more than one form of social organization was observed in the same population. Populations where 2 or more forms of group-living but no other forms of social organization occurred were categorised as “variable group”. Otherwise, to reduce the number of IVSO categories for the statistical analysis, we combined the four forms of group living as “group” when the population had both group-living and non-group-living social units. Therefore, IVSO consisted of five categories: (1) solitary + pair (SP), (2) solitary + group (SG), (3) pair + group (PG), (4) solitary + pair + group (SPG), (5) variable group (VG). Together with the 3 non-IVSO categories: solitary (S), pair (P) and stable group (G, only one form of group-living reported), we have 8 combinations of social organization, but as SG did not occur in any population, this was effectively reduced to 7. If males and females live separately in two different forms of social organization (for example, group-living females and solitary males), this population was classified as having a sex-specific social organization and not IVSO, since there was no variation within either sex. Fission-fusion groups characterized by temporal variation in group size and composition [2] are common in some kangaroo species [3]. All observed variation in social organization in a fission-fusion population was recorded. Environmental disruption events, such as the accidental death of a group member, represent external incidents that can change the social organization of the social unit. Three cases of social organization change due to environmental disruptors were not considered for further analysis because they do not represent a change in social organization in response to ancestrally relevant conditions.

To better estimate the ancestral state and the possible ecological factors that may have an influence on the evolution of social organization, we obtained body mass from the Handbook Mammals of the World [4] and centred it to the body mass of Australia’s oldest known marsupial fossil [5]. By comparing the M2 mesiodistal length, we estimated the body mass of the 30 million years old ancestral species Djarthia murgonensis would be similar to Antechinus stuartii as 37.75g, which is smaller than most of the modern species in our database and much smaller than the mean (4871g). The number of studies per population was recorded to control for research effort. We determined the habitat type(s) in which the study took place, categorized and recorded based on IUCN classification (www.iucn.org) as desert, forest, rocky areas, savannah, grassland, shrubland, wetlands or artificial. Habitat heterogeneity was then determined as the maximum number of habitats occupied per population.

The ancestral state of social organization and climate in Australian marsupials: To test how climate would affect social organization, we focused on Australian marsupial species to control for other environmental variables that differ between Australia and South America, such as competition with placentals that mainly occurs in South America but not Australia (apart from small rodents and bats). We obtained long-term climate data from the online dataset of the Australian Bureau of Meteorology (http://www.bom.gov.au/climate). For each studied population, we obtained local climate data at the GPS coordinates reported in the articles. Climate data were obtained for 51 Australian marsupial species.

High-resolution (0.05x0.05 degree) grids downloaded from the climate dataset were converted to raster grids in QGIS 3.10. Based on the monthly precipitation and monthly mean maximum temperature data from 1910 to 2019, we calculated six variables to represent climate conditions: annual mean precipitation (mm), annual mean maximum temperature (°C) and coefficient of variation to represent within-year variation (seasonality) and between-year variation (predictability) for both precipitation and temperature (for details see Supplementary materials S1). Six maps were generated in QGIS, one for each climate variable. For populations with precise GPS locations, data were directly obtained from the climate maps. When the location was not specified in the paper (one population of Distoechurus pennatus in the state of Victoria), we ran “zonal statistics” analysis to obtain the mean value of the area.

Statistical analysis: Phylogenetic comparative analyses were conducted by R v.3.6.1, using the R packages brms [6, 7], RStan [8] and Rethinking [9]. All R codes and data are available at https://github.com/JingyuQ/MarsupialSO.

Climate data are often correlated with each other. Thus, we first performed a principal component analysis (PCA) to reduce the six climate variables (see above) to a smaller number of components. The first two principal components (PC1 and PC2) explained 82% of the variation. PC 1 was positively related with annual temperature (eigenvector=0.506) and within-year variation of precipitation (eigenvector =0.449), and negatively related with within-year variation of temperature (eigenvector=-0.497) and between-year variation of temperature (eigenvector=-0.494). High PC1 values match the tropical rainforest climate in Australia. PC2 was positively related to annual precipitation (eigenvector=0.693) and negatively with the between-year variation of precipitation (eigenvector=-0.639). Low PC2 values match the desert climate in central Australia.

We used Bayesian generalized linear mixed-effects models (GLMMs) to control for phylogeny and estimate the associations between social organization and predictor variables [10]. The probability of each kind of social organization was modelled using a categorical likelihood, allowing gradual changes in the probability of each kind of social organization along the phylogeny; the intercept of such a model represents the phylogenetically-controlled mean of extant species, and, in the absence of any directional trends, the ancestral state [10, 11]. This approach corresponds to a polygenic model of trait inheritance, as opposed to alternative approaches to inferring ancestral states of categorical traits based on nucleotide substitution models [12, 13]; in substitution models, evolutionary change is not gradual but occurs in “jumps” between states, which has been deemed less plausible than polygenic inheritance [10, 14] (and in our opinion rightfully so). Furthermore, these models often have many more free parameters, require treating the species rather than the population as the unit of observation (thus failing to account for intra-specific variation), and cannot include covariates when inferring ancestral states. We therefore prefer the quantitative genetic approach implemented by our GLMMs. This does not mean that social organization is directly genetically inherited, but simply that whatever individual-level traits contribute to social organization follow a polygenic rather than single-gene pattern of inheritance (which is widely accepted to be the case for social and behavioural traits, see e.g. [15]).

The phylogenetic history and its uncertainty were represented by a sample of 100 phylogenetic trees, downloaded from VertLife (http://vertlife.org/phylosubsets/) [16]. We created two models, the first one was to estimate the ancestral state of modern marsupials, therefore included data from all studied marsupial species: social organization (with 7 categories) ~ species intercept + covariance for phylogeny + habitat heterogeneity + number of studies + body mass. The second model aimed to estimate ecological (climate and habitat) effects on Australian marsupial’s social organization: social organization ~ species intercept + covariance for phylogeny + habitat heterogeneity + number of studies + body mass + climate PC 1 + climate PC 2. Due to the limitations of climate data, this model only considered Australian marsupials. Both models were run at the population level and included phylogeny and species as random factors (as indicated above). The number of studies per population was considered a predictor of the occurrence of IVSO. Phylogenetic signal (λ) was calculated as the proportion of variance captured by the phylogenetic random effect [47], representing the tendency of related species to resemble each other more than species drawn at random from the same tree [48]. For analysis, the social organization of populations showing only one form of group living was categorised as a stable group. For more details on model structure, see the PDF “Model details”.

1.     Makuya L, Olivier C-A, Schradin C. In press. Field studies need to report essential information on social organization – independent of the study focus. Ethology n/a. (doi:10.1111/eth.13249)

2.     Aureli F et al. 2008 Fission‐Fusion Dynamics: New Research Frameworks. Current Anthropology 49, 627–654. (doi:10.1086/586708)

3.     Best EC, Seddon JM, Dwyer RG, Goldizen AW. 2013 Social preference influences female community structure in a population of wild eastern grey kangaroos. Animal Behaviour 86, 1031–1040. (doi:10.1016/j.anbehav.2013.09.008)

4.     Wilson, D,E, & Mittermeier, R,A. 2015 Handbook of the Mammals of the World Volumne 5: Monotremes and Marsupials. Lynx Edicions, Barcelona.

5.     Beck RMD, Godthelp H, Weisbecker V, Archer M, Hand SJ. 2008 Australia’s Oldest Marsupial Fossils and their Biogeographical Implications. PLoS ONE 3, e1858. (doi:10.1371/journal.pone.0001858)  

6.     Bürkner P. C. 2017 brms: An R Package for Bayesian Multilevel Models Using Stan. Journal of Statistical Software 80, 1–28. (doi:10.18637/jss.v080.i01)  

7.      Bürkner P. C. 2018. Advanced Bayesian Multilevel Modeling with the R Package brms. The R Journal 10(1), 395-411. (doi.org/10.32614/RJ-2018-017)  

8.     Stan Development Team. 2020 RStan: the R interface to Stan. R package version 2.21.2. (http://mc-stan.org/)  

9.     McElreath R. 2020 Statistical rethinking: A Bayesian course with examples in R and Stan. Chapman and Hall/CRC.

10.  Hadfield JD, Nakagawa S. 2010 General quantitative genetic methods for comparative biology: phylogenies, taxonomies and multi-trait models for continuous and categorical characters. Journal of Evolutionary Biology 23, 494–508. (doi:10.1111/j.1420-9101.2009.01915.x)  

11.  Lynch M. 1991 Methods for the Analysis of Comparative Data in Evolutionary Biology. Evolution 45, 1065–1080. (doi:10.1111/j.1558-5646.1991.tb04375.x)  

12.  Pagel M. 1994 Detecting correlated evolution on phylogenies: a general method for the comparative analysis of discrete characters. Proceedings of the Royal Society of London. Series B: Biological Sciences 255, 37–45. (doi:10.1098/rspb.1994.0006)  

13.  Pagel M, Meade A. 2006 Bayesian Analysis of Correlated Evolution of Discrete Characters by Reversible‐Jump Markov Chain Monte Carlo. The American Naturalist 167, 808–825. (doi:10.1086/503444)  

14.  Felsenstein J. 2005 Using the quantitative genetic threshold model for inferences between and within species. Philosophical Transactions of the Royal Society B: Biological Sciences 360, 1427–1434. (doi:10.1098/rstb.2005.1669)  

15.  Chabris CF, Lee JJ, Benjamin DJ, Beauchamp JP, Glaeser EL, Borst G, Pinker S, Laibson DI. 2013 Why It Is Hard to Find Genes Associated With Social Science Traits: Theoretical and Empirical Considerations. Am J Public Health 103, S152–S166. (doi:10.2105/AJPH.2013.301327)

16.  Inferring the mammal tree: species-level sets of phylogenies for questions in ecology, evolution, and conservation; PLOS Biology, 2019; Upham NS, Esselstyn JA, Jetz W., doi:10.5061/dryad.tb03d03

17.  Nakagawa S, Schielzeth H. 2012 The mean strikes back: mean–variance relationships and heteroscedasticity. Trends in Ecology and Evolution  

18.  Münkemüller T, Lavergne S, Bzeznik B, Dray S, Jombart T, Schiffers K, Thuiller W. 2012 How to measure and test phylogenetic signal. Methods in Ecology and Evolution 3, 743–756. (doi:10.1111/j.2041-210X.2012.00196.x)