Investigating COVID-19 transmission and mortality differences between indigenous and non-indigenous populations in Mexico.

OBJECTIVES: Indigenous populations have been disproportionately affected during pandemics. We investigated COVID-19 mortality estimates among indigenous and non-indigenous populations at national and sub-national levels in Mexico.
METHODS: We obtained data from the Ministry of Health, Mexico, on 2,173,036 laboratory-confirmed RT-PCR positive COVID-19 cases and 238,803 deaths. We estimated mortality per 1000 person-weeks, mortality rate ratio (RR) among indigenous vs. non-indigenous groups, and hazard ratio (HR) for COVID-19 deaths across four waves of the pandemic, from February 2020 to March 2022. We also assessed differences in the reproduction number (Rt).
RESULTS: The mortality rate among indigenous populations of Mexico was 68% higher than that of non-indigenous groups. Out of 32 federal entities, 23 exhibited higher mortality rates among indigenous groups (P < 0.05 in 13 entities). The fourth wave showed the highest RR (2.40). The crude HR was 1.67 (95% CI: 1.62, 1.72), which decreased to 1.08 (95% CI: 1.04, 1.11) after controlling for other covariates. During the intense fourth wave, the Rt among the two groups was comparable.
CONCLUSION: Indigenous status is a significant risk factor for COVID-19 mortality in Mexico. Our findings may reflect disparities in non-pharmaceutical (e.g., handwashing and using facemasks), and COVID-19 vaccination interventions among indigenous and non-indigenous populations in Mexico.

Introduction

Globally, indigenous populations and ethnic minorities tend to suffer worse health outcomes than non-indigenous populations. For example, from 2010-2012 the life expectancy of indigenous Australians (Aboriginal and Torres Strait Islander) (69.1 years for males and 73.7 years for females) was around 10 years lower than that of non-indigenous Australians (79.7 years for males and 83.1 years for females) (Australian Institute of Health Welfare, 2014). A variety of factors such as transgenerational adverse effects of colonization, racism, lower socio-economic status, and lower levels of education, have contributed to existing health disparities among the indigenous population (Durey ; Hajizadeh ; Nazroo, 2003). Moreover, barriers within the healthcare setting, including lack of access to culturally appropriate healthcare services, poor health literacy, distance to medical centers, communication problems, and low health insurance coverage also increase negative health trajectories (Daws ; Walsh and Kangaharan, 2017).

Studies from many countries have reported disproportionately higher infection rates, hospitalizations, and deaths among indigenous and ethnic-minority groups from the beginning of the COVID-19 pandemic (Mallard ; Power ; Sharma and Bhaskar, 2021; Wiemers ; Yashadhana ). However, a few countries have reported the opposite early in the COVID-19 pandemic. For example, a review article from November 2020, covering COVID-19 studies reported that in six out of nine countries (including Mexico) that publicly reported mortality data among indigenous peoples, the rate of infection per 100,000 population was lower among indigenous groups (except Brazil, Peru, and USA), and the case fatality ratio was lower in three countries (Mallard ). However, it was probably too early in the pandemic to assess this comparison since some indigenous communities tend to be more isolated (Tripp, 2022). Moreover, there are major gaps in reporting data on COVID-19 infection and outcomes for indigenous populations and many countries do not have publicly available data with socio-demographic profiles, leading to a dearth of studies that have empirically examined COVID-19 in indigenous populations vs. non-indigenous populations (Alves ). Therefore, there is a dire need to conduct studies that evaluate COVID-19 outcomes among indigenous and non-indigenous populations.

The COVID-19 pandemic widened existing health disparities among these groups as a result of the higher prevalence of underlying comorbid conditions and multimorbidity compounded by poor access to health care services, and unequal impact of lockdown measures (e.g., higher unemployment, higher rates of mental health problems, more loss of schooling) (Katikireddi ; Kirby, 2020; Yashadhana ). Prior work indicates that indigenous status may influence pandemic outcomes through three key mechanisms: a) increased exposure (at work, at home, possibly through multigenerational living), b) increased medical susceptibility due to higher prevalence of non-communicable diseases, and c) limited access to care and or health literacy resulting in social disparities (Quinn and Kumar, 2014).

Mexico is a country in the Americas with the largest indigenous population (IWGIA, 2021). Mexico's indigenous population is highly vulnerable to the effects of the COVID-19 pandemic due to factors such as marginalization, discrimination, violence, land dispossession, and poor access to health services, social security, education and adequate housing (Díaz de León-Martínez ; IWGIA, 2021; Consejo Nacional de Evaluación de la Politica de Desrrollo Social, 2019). However, the COVID-19 prevention and mitigation strategies implemented in the country for indigenous populations are the same as those for the general population (Díaz de León-Martínez ). Preliminary analyses from Mexico have started to shed light on the severe health impacts of COVID-19 among indigenous populations compared to non-indigenous counterparts (Argoty-Pantoja ; Ibarra-Nava ; Mallard ). In this study, we aimed to assess the mortality impact of COVID-19 among indigenous populations by quantifying the mortality rate at sub-national level and by comparing the rates between indigenous and non-indigenous populations across the waves of the COVID-19 pandemic. In addition, we compared the hazard ratios for death between indigenous and non-indigenous populations controlling for socio-demographic, comorbidity-related, and health service utilization-related factors. Furthermore, we estimate the reproduction number (R) among the two groups. These detailed analyses can help guide intervention strategies that effectively address the impact of infectious disease emergencies such as the COVID-19 pandemic among vulnerable populations.

Methods

Study setting

Mexico is inhabited by 128.97 million people (projected as of 2021) and is divided into 31 states and the federal district (Mexico City) collectively referred to as federal entities (CONAPO, 2019). According to the intercensal survey of 2015 in Mexico, 21.50% of the total population self-identified as indigenous, 65% of whom were concentrated in 8 of the 32 entities (National Institute of Indigenous Peoples, 2017) as shown in table 2.

Data source

The dataset for this study was retrieved from the website (Secretaría de Salud, 2022) of the General Directorate of Epidemiology of the Ministry of Health of Mexico which maintains an open-source repository of patients classified as ‘suspected cases of viral respiratory disease’ identified at medical facilities in Mexico through the Viral Respiratory Diseases Epidemiological Surveillance System (Secretaría de Salud, 2020). From this dataset, we obtained data on the selected variables: age (EDAD), sex (SEXO), state of patient's residence (ENTIDAD_RES), type of care the patient received (TIPO_PACIENTE), date of symptom onset (FECHA_SINTOMAS), date of admission to care (FECHA_INGRESO), date of death (FECHA_DEF), lab result (RESULTADO_LAB), patient's self-identification as indigenous (INDIGENA), and variables on the diagnosis of different conditions such as diabetes, hypertension, and obesity which have been specified in Table 1 . We used the variable ‘lab result’ to select the confirmed cases of COVID-19 for our study. We categorized five groups based on age and created a composite variable ‘comorbidity category’ based on the presence or absence of nine different comorbid conditions (Table 1).

Table 1

Descriptive analysis of socio-demographic, comorbidity, and case management-related characteristics of COVID-19 positive cases among indigenous and non-indigenous populations.

	SARS-Cov 2 positive (n = 2,173,036)		Mortality among SARS-Cov2 positive (n = 238,803)		Proportion of deaths among cases
	Indigenous population (n = 21,896)	Non-indigenous population (n = 2,151,140)	Indigenous population (n = 3831)	Non-indigenous population (n = 234,972)	Indigenous	Non-indigenous	P-valuea
Sex n (%)
Female	10,358 (47.31)	1,079,312 (50.17)	1473 (38.45)	89,604 (38.13)	14.22	8.30	<0.0001
Male	11,538 (52.69)	1,071,828 (49.83)	2358 (61.55)	145,368 (61.87)	20.43	13.56	<0.0001
Age group n (%)
Less than 18 years	861 (3.93)	90,048 (4.19)	26 (0.68)	896 (0.38)	3.02	0.99	<0.0001
18-44 years	9113 (41.62)	107,3787 (49.92)	335 (8.74)	24,214 (10.31)	3.68	2.25	<0.0001
45-54 years	3914 (17.88)	407,633 (18.95)	555 (14.49)	36,798 (15.66)	14.18	9.03	<0.0001
55-64 years	3453 (15.77)	290,059 (13.48)	945 (24.67)	58,293 (24.81)	27.37	20.09	<0.0001
65 and above years	4555 (20.80)	289,613 (13.46)	1970 (51.42)	114,771 (48.84)	43.25	39.63	<0.0001
Type of care n (%)
Outpatient	14,179 (64.76)	1,650,958 (76.75)	254 (6.63)	15,733 (6.70)	1.79	0.95	<0.0001
Hospitalized	7717 (35.24)	500,182 (23.25)	3577 (93.37)	219,239 (93.30)	46.35	43.83	<0.0001
Patient required admission to intensive care unit n (%)	732 (3.34)	41,105 (1.91)	431 (11.25)	23,809 (10.13)	58.88	57.92	0.603
Comorbidity
Pneumonia, n (%)	6011 (27.45)	364,845 (16.96)	3023 (78.91)	171,513 (72.99)	50.29	47.01	<0.0001
Diabetes, n (%)	4156 (18.98)	300,200 (13.96)	1425 (37.20)	86,625 (36.87)	34.27	28.86	<0.0001
COPD, n (%)	625 (2.85)	26,266 (1.22)	266 (6.94)	10,240 (4.36)	42.56	38.99	0.070
Asthma, n (%)	521 (2.38)	48,438 (2.25)	96 (2.51)	4231 (1.79)	18.43	8.73	<0.0001
Hypertension, n (%)	4521 (20.65)	387,275 (18.00)	1520 (39.68)	104,644 (44.53)	33.62	27.02	<0.0001
Other disease, n (%)	465 (2.12)	44,684 (2.08)	146 (3.81)	11,577 (4.93)	31.40	25.91	0.007
Cardiovascular disease, n (%)	438 (2.00)	36,111 (1.68)	167 (4.36)	11,859 (5.05)	38.13	32.84	0.019
Obesity, n (%)	3910 (17.86)	311,392 (14.48)	907 (23.68)	51,063 (21.73)	23.20	16.40	<0.0001
Chronic kidney failure, n (%)	417 (1.90)	37,271 (1.73)	197 (5.14)	16,097 (6.85)	47.24	43.19	0.097
Comorbidity category n (%)
No comorbidity	9877 (45.11)	1,239,373 (57.61)	245 (6.40)	20,485 (8.72)	2.48	1.65	<0.0001
1-2 comorbidities	9658 (44.11)	743,630 (34.57)	2369 (61.84)	138,669 (59.02)	24.53	18.65	<0.0001
Three or more comorbidities	2361 (10.78)	168,137 (7.82)	1217 (31.77)	75,818 (32.27)	51.55	45.09	<0.0001
Time related variablesc					P-valueb
Time from symptom onset to seeking care (in days)d(n=2,169,701)	4.41 (3.20)	4.07 (3.17)	-	-	<0.0001		<0.0001
Time from symptom onset to death (in days)e(n=237,401)	-	-	13.25 (8.22)	13.89 (8.53)	<0.0001		<0.0001
Time from seeking care to death (in days)f(n=236,134)	-	-	7.63 (7.08)	8.45 (7.34)	<0.0001		<0.0001

Two-sample test for equality of proportions of mortality out of total cases among indigenous and non-indigenous groups for the given characteristic

Non-parametric Wilcoxon two-sample test for difference in mean time among indigenous and non-indigenous

Mean (SD)

Includes cases with time from 0-18 days

Includes cases with time from 0-52 days

Includes cases with time from 0-39 days

Definition of confirmed case of COVID-19

Laboratory-confirmed cases were defined as a positive RT-PCR test (Ibarra-Nava ).

Definition of indigenous populations

The use of self-identification as a proxy for ethnic classification has been used by the Mexican census since 2000 and this method is in line with the spirit of international legislation that considers the ability of the indigenous population to identify their ethnicity as a fundamental right (ILO, 1989). In this study, we used the variable ‘INDIGENA’ from the dataset to measure indigenous status based on patient self-identification. Records, where the indigenous status was missing, were omitted from further analysis.

Figure 1 displays the strategy used to determine inclusion of cases in our study. After excluding cases without information about indigenous status, a total of 2,173,036 COVID-19 positive cases, and 238,803 COVID-19 deaths remained.

Figure 1

Flowchart showing the inclusion of cases and deaths in the study.

Data analysis

The earliest date of symptom onset for a COVID-19 positive case in our dataset was 19 February 2020. Therefore, we performed the statistical analysis using the surveillance data from 19 February 2020, to 25 March 2022.

a. Descriptive analysis of COVID-19 positive cases and deaths: To test the null hypothesis of no difference in the proportion of COVID-19 deaths among indigenous and non-indigenous positive cases for socio-demographic, comorbidity, and case management-related characteristics, we performed two-sample tests for equality of binomial proportions at 0.05 level of significance. Similarly, to test the null hypothesis of no difference in the average days from ‘symptom onset to seeking care’, ‘seeking care to death’ and ‘symptom onset to death’ between indigenous and non-indigenous populations, we used the Wilcoxon two-sample test for difference in medians. For the analysis of these time-related variables, we excluded cases with negative days, and extreme outliers falling above three times the interquartile range.

b. Estimation of instantaneous reproduction number: R is the expected number of secondary infections occurring at time t, divided by the number of infected individuals, each scaled by their relative infectiousness at time t. We used EpiEstim R Package (Cori ) to estimate R from the curves of daily incidence of COVID-19 cases by date of symptom onset for indigenous and non-indigenous populations separately. We obtained the average R estimates over a weekly time interval for the entire study period by using the 7-day sliding window method. We also used the non-overlapping time window method to estimate the average R for the early ascending phase for each of the COVID-19 waves. For the R estimation, the distribution of serial interval was parametrically defined using a mean of 4.6 days and SD of 5.55 days (Ofori ; You ). We report the mean and 95% credible interval (CrL).

c. Estimation of person-time mortality rate: The person-time mortality rate was defined as the ratio of the number of deaths among COVID-19 cases and the person-time at risk of death during the study period. The mortality rate was estimated separately for indigenous and non-indigenous populations. The person-time at risk was expressed per 1000 person-weeks based on the date from symptom onset to date of death (Argoty-Pantoja ). For non-deaths, person-time at risk was the time between date of symptom onset until the last date of the study period. We also estimated the person-time mortality rate for each of the 32 federal entities. The 95% confidence intervals (CIs) and P-values were estimated using open-source statistics for public health (Sullivan ). The P-values were based on z-score tests and the 95% CIs were based on the Taylor series method. We also estimated the mortality rate for four different waves of COVID-19 as well as for the years 2020 and 2021 in Mexico. For the wave-wise comparison, person-time for non-deaths was estimated using the time elapsed between symptom onset date and the end date of the given wave. We also reported the ratio of COVID-19 mortality rate per 1000 person-week among indigenous populations vs. non-indigenous populations (rate ratio = RR).

We defined the epidemiological waves based on the descriptive analysis of time series of COVID-19 cases for Mexico according to the standard definition of epidemiological weeks used by the CDC (CMMCP, 2022) as follows:

First wave- 19 February 2020 to 3 October 2020; second wave- 4 October 2020 to 29 May 2021; third wave- 30 May 2021 to 18 December 2021; fourth wave- 19 December 2021 to 25 March 2022.

d. Estimation of hazard ratio (HR) for COVID-19 death and indigenous status: Multivariate Cox proportional hazards regression models were used for these analyses. The general equation of the model iswhere  is the hazard function, is the baseline hazard, X is a variable indicating indigenous status, and is a vector of model covariates including age, sex, number of comorbidities, type of care received (outpatient vs. hospitalization), and days from symptom onset to seeking care. We used log-log survival curves to assess the constant proportionality assumption of the Cox hazard model. Because the curves showed parallel lines for indigenous and non-indigenous populations, the constant proportionality assumption was reasonable.

Results

Figure 2 shows the daily time-series of COVID-19 cases from 19 February 2021 to 25 March 2022 (panels A and D), and R among indigenous and non-indigenous populations in Mexico (panels B, C, E and F). The mean R was highest among both the groups during the early ascending phase of the fourth wave (indigenous: 1.426 (95% CrL: 1.405, 1.447), non-indigenous: 1.535 (95% CrL: 1.521, 1.549)) (Supplementary Figure 2).

Figure 2

(A.) Daily incidence of COVID-19 cases among non-indigenous from 19 February 2020 - 25 March 2022 (B.) Time-varying R of COVID-19 cases among non-indigenous populations for the study period (C.) Mean R for the early ascending phase of four waves among non-indigenous (30 days from 5 April 2020 - 4 May 2020 for first wave, 30 days from 29 September 2020 – 28 October 2020 for second wave, 30 days from 27 June 2021 - 26 July 2021 for the third wave and 15 days from 23 December 2021 - 6 January 2022 for the fourth wave) (D.) Daily incidence of COVID-19 cases among Indigenous from 19 February 2020 -25 March 2022 (E.) Time-varying R of COVID-19 cases among indigenous populations for the study period (F.) Mean R for the early ascending phase of four waves among indigenous (30 days from 15 April 2020 - 14 May 2020 for first wave, 30 days from 15 November 2020 - 14 December 2020 for second wave, 30 days from 11 July 2021 - 9 August 2021 for the third wave and 15 days from 27 December 2021 - 10 January 2022 for the fourth wave). Dates indicate date of symptom onset. The shaded gray areas in B, C, E and F indicate 95% CrL.

CrL, credible interval; R, reproduction number

Table 1 summarizes socio-demographic, comorbidity, and case management-related characteristics among COVID-19 cases and deaths for the indigenous and non-indigenous populations. The proportion of deaths among the indigenous group was higher than that for non-indigenous for both males and females, across age groups, both hospitalized and outpatient cases, as well as the great majority of comorbidities (P <0.05). Compared to non-indigenous populations, the average time from symptom onset to seeking care was significantly higher for the indigenous populations (4.41 vs. 4.07 days; P<0.01) whereas both the average time from symptom onset to death (13.25 vs. 13.89 days; P<0.01) and the average time from seeking care to death (7.63 vs. 8.45 days; P<0.01) were significantly lower for the indigenous populations (Table 1). These findings indicate longer delays in care seeking for the indigenous population leading to lower average survival time compared to the non-indigenous populations.

Table 2 shows the federal entity-level COVID-19 mortality rate per 1000 person-weeks among indigenous and non-indigenous populations. At the national level, the mortality rates among indigenous and non-indigenous populations were significantly different (3.25 vs. 1.94 per 1000 person-weeks; P-value<0.05). Of the 32 federal entities, 23 exhibited higher COVID-19 mortality rate among indigenous compared to non-indigenous groups, with this difference being statistically significant in 13 entities. Figure 3, Figure 4 show the mortality rates by indigenous status and the RR across federal entities, respectively. We also assessed the association between the proportion of indigenous populations and the RR at the entities level (Pearson correlation coefficient = 0.37; P-value = 0.036; 95% CI: 0.02, 0.64).

Table 2

COVID-19 mortality rate per 1000 person-weeks among indigenous and non-indigenous population, Mexico, 19 February 2020 to 25 March 2022.

Region/entities	Proportion of self-identified indigenous population	Total deaths		Mortality rate (95% CI)
		Indigenous	Non-indigenous	Indigenous	Non-indigenous	P-value
National	21.49	3818	233,750	3.25 (3.15, 3.35)	1.94 (1.93, 1.94)	<0.0001
Aguacalientes	0.60	7	2544	3.19 (1.28, 6.56)	1.69 (1.63, 1.76)	0.089
Baja California	1.10	83	8862	3.96 (3.15, 4.90)	3.24 (3.17, 3.31)	0.069
Baja California Sur	0.40	10	1881	2.62 (1.25, 4.82)	0.97 (0.93, 1.02)	0.001
Campeche	1.56	79	1936	4.78 (3.78, 5.95)	2.51 (2.40, 2.63)	<0.0001
Coahuila	0.80	12	7143	1.68 (0.87, 2.93)	1.69 (1.66, 1.73)	0.976
Colima	0.57	6	1725	2.66 (0.97, 5.79)	2.53 (2.41, 2.65)	0.902
Chiapas	7.34	34	2008	2.00 (1.38, 2.79)	2.48 (2.37, 2.59)	0.206
Chihuahua	1.56	71	6595	3.12 (2.44, 3.93)	2.37 (2.31, 2.42)	0.020
Mexico City	3.05	127	26,972	1.44 (1.20, 1.72)	1.28 (1.27, 1.30)	0.183
Durango	0.54	12	2927	1.26 (0.65, 2.19)	1.29 (1.25, 1.34)	0.915
Guanajuato	2.08	11	12,344	1.01 (0.50, 1.81)	1.50 (1.47, 1.53)	0.188
Guerrero	4.66	144	4400	2.21 (1.86, 2.60)	1.95 (1.90, 2.01)	0.144
Hidalgo	4.03	315	6810	3.80 (3.40, 4.25)	3.27 (3.19, 3.35)	0.008
Jalisco	3.40	65	14,691	3.61 (2.78, 4.60)	2.94 (2.90, 2.99)	0.102
Mexico	10.71	224	31,415	3.36 (2.93, 3.82)	2.93 (2.90, 2.96)	0.042
Michoacan	4.94	190	6249	3.58 (3.08, 4.12)	2.17 (2.11, 2.22)	<0.0001
Morelos	2.08	30	3456	4.86 (3.27, 6.93)	3.83 (3.70, 3.96)	0.193
Nayarit	1.02	31	1899	2.91 (1.98, 4.13)	2.89 (2.76, 3.02)	0.966
Nuevo Leon	1.37	15	13,063	1.46 (0.82, 2.41)	1.71 (1.68, 1.74)	0.549
Oaxaca	10.15	436	3476	4.02 (3.65, 4.41)	1.27 (1.23, 1.31)	<0.0001
Puebla	8.47	313	11,266	7.51 (6.70, 8.39)	2.41 (2.37, 2.46)	<0.0001
Queretaro	1.52	13	4151	1.20 (0.64, 2.05)	1.19 (1.15, 1.23)	0.981
Quintana Roo	2.60	205	3452	5.44 (4.72, 6.24)	2.63 (2.54, 2.72)	<0.0001
San Luis Potosi	2.45	139	5188	1.34 (1.13, 1.58)	1.38 (1.34, 1.41)	0.768
Sinaloa	1.48	27	7078	3.47 (2.29, 5.06)	3.05 (2.97, 3.12)	0.494
Sonora	1.98	90	6795	4.81 (3.87, 5.91)	1.87 (1.83, 1.92)	<0.0001
Tabasco	2.40	76	5379	0.99 (0.78, 1.24)	0.92 (0.89, 0.94)	0.528
Tamaulipas	0.84	10	6125	1.35 (0.64, 2.48)	1.43 (1.39, 1.46)	0.854
Tlaxcala	1.25	18	2721	3.12 (1.85, 4.92)	2.33 (2.24, 2.42)	0.217
Veracruz	9.24	190	12,879	5.78 (4.99, 6.66)	3.19 (3.13, 3.24)	<0.0001
Yucatan	5.34	822	4935	4.04 (3.77, 4.33)	2.17 (2.11, 2.23)	<0.0001
Zacatecas	0.47	13	3385	2.06 (1.10, 3.52)	1.67 (1.62, 1.73)	0.451

Figure 3

COVID-19 mortality rate per 1000 person-weeks in indigenous and non-indigenous populations in Mexico, February 2020 to March 2022.

Figure 4

Rate ratio (RR) of COVID-19 mortality per 1000 person-weeks among indigenous and non-indigenous populations by federal entities in Mexico.

A stratified analysis of the mortality rate and the RR across four pandemic waves is given in Table 3 . The RR was lowest in the first wave (1.55) and highest in the fourth wave (2.40). Similarly, the RR was higher in 2021 compared to 2020.

Table 3

Estimation of COVID-19 mortality rate per 1000 person-week among indigenous and non-indigenous groups across four pandemic waves and by year in Mexico.

Wave	Deaths indigenous	Deaths non-indigenous	Person-weeks indigenous	Person-weeks non-indigenous	Mortality rate per 1000 person-weeks indigenous (R1)	Mortality rate per 1000 person-weeks non-indigenous (R2)	Rate ratio (R1/R2)
First wave19 Feb 2020 – 3 Oct 2020	1679	83,002	99,912.42	7,693,327	16.71	10.79	1.55
Second wave4 Oct 2020 – 29 May 2021	985	95,517	107,946.1	17,323,003	9.12	5.51	1.65
Third wave30 May 2021 – 18 Dec 2021(First case of Delta variant reported in 15 July 2021)	926	38,627	68,398.43	5,309,029	13.54	7.27	1.86
Fourth wave19 Dec 2021 – 25 March 2022(First case of Omicron variant reported on 3 Dec 2021)	141	9666	10,932.43	1,800,412	12.90	5.37	2.40
2020(19 Feb 2020 – 31 Dec 2021)	2055	127,131	902,822.9	93,530,292.14	2.28	1.36	1.67
2021(1 Jan 2021 – 31 Dec 2021)	1533	85,195	262,511.9	25,549,998.29	5.84	3.33	1.75
Overall (19 Feb 2020 – 25 March 2022)	3818	233,750	1,175,598	120,703,675.1	3.25	1.94	1.68

Note: The overall person-time was calculated as the difference between the date of symptom onset and death for events (death). For non-death cases (censored) the person-time was calculated as the difference between date of symptom onset and cutoff date of 25 March 2021. For wave-wise estimation, person-time was calculated for the interval period only. Person-week was calculated by dividing total person-days by 7.

Table 4 shows the results from fitting different survival regression models. Based on the goodness of fit of the models as deemed by their Akaike information criterion (AIC) values, we selected model 6 as our final model. In the unadjusted model, the hazard of death among indigenous populations was 67.1% (HR = 1.67; 95% CI: 1.62, 1.72) higher than that for non-indigenous populations. There was a statistically significant interaction between the type of care and the time from symptom onset to seeking care (care-seeking delay). After adjusting for sex, age, the number of comorbidities, and the interaction between type of care and care-seeking delay, the HR for indigenous groups compared to non-indigenous decreased to 1.08 and remained statistically significant. We used the variables in model 6 to fit survival regression models for each of the waves separately. Table 5 shows the results of these models. Tables with the estimated beta values for model 6 for the overall period, waves, and years are given in the supplementary file.

Table 4

Results of fitting a taxonomy of multiple survival regression models, February 2020 to March 2022.

	Hazard ratio (95% CI)
	Model 1	Model 2	Model 3	Model 4	Model 5	Model 6
Indigenous(Comparison = Non-indigenous)	1.67 ***(1.62, 1.72)	1.24***(1.20, 1.28)	1.61***(1.56, 1.67)	1.12***(1.09, 1.16)	1.08***(1.05, 1.11)	1.08***(1.04, 1.11)
Sex(Comparison = Female)		1.55***(1.54, 1.56)			1.27***(1.26, 1.28)	1.27***(1.26, 1.28)
Age		1.07***(1.07, 1.07)			1.03***(1.03, 1.03)	1.03***(1.03, 1.03)
Time in days from symptom onset to seeking care			1.10***(1.10, 1.10)	0.99***(0.99, 0.99)	0.98***(0.98, 0.98)	-
Type of care(Comparison = Outpatient)				59.67***(58.71, 60.65)	24.55***(24.13, 24.99)	-
1-2 comorbidity(Comparison = No comorbidity)					2.30***(2.26, 2.34)	2.28***(2.25, 2.32)
3 or more comorbidity(Comparison = No comorbidity)					2.96***(2.91, 3.01)	2.94***(2.90, 2.99)
(Type of care) * (time in days from symptom onset to seeking care)
Time in days from symptom onset to seeking care (for outpatients)						1.07***(1.06, 1.07)
Time in days from symptom onset to seeking care (for hospitalized)						0.98***(0.97, 0.98)
Type of care, when time from symptom onset to care=0(Comparison = Outpatient)						36.29***(35.28, 37.33)
Akaike information criterion	6,879,429.7	6,552,629.9	6,849,913.9	6,260,272.2	6,169,794.7	6,168,502.6

*P<0.05, **P<0.01, ***P<0.001

CI, confidence interval

Table 5

Results of multiple survival regression models across COVID-19 waves.

	Hazard ratio (95% CI)
	First wave	Second wave	Third wave	Fourth wave
Indigenous(Comparison = Non-indigenous)	1.09***(1.03, 1.14)	1.04(0.98, 1.11)	1.13**(1.06, 1.21)	1.03(0.87, 1.22)
Sex(Comparison = Female)	1.32***(1.30, 1.34)	1.26***(1.24, 1.28)	1.19***(1.17, 1.22)	1.25***(1.20, 1.30)
Age	1.03***(1.03, 1.03)	1.03***(1.03, 1.03)	1.02***(1.02, 1.02)	1.03***(1.03, 1.03)
1-2 comorbidity(Comparison = No comorbidity)	2.69***(2.62, 2.77)	2.06***(2.01, 2.11)	1.93***(1.87, 2.00)	1.88***(1.76, 2.00)
3 or more comorbidity(Comparison = No comorbidity)	3.69***(3.58, 3.80)	2.58***(2.51, 2.65)	2.41***(2.32, 2.50)	2.25***(2.09, 2.41)
(Type of care) * (time in days from symptom onset to seeking care)
Time in days from symptom onset to seeking care (for outpatients)	1.04***(1.03, 1.04)	1.07***(1.07, 1.08)	1.03**(1.01, 1.05)	1.11***(1.08, 1.15)
Time in days from symptom onset to seeking care (for hospitalized)	0.96***(0.96, 0.96)	0.97***(0.97, 0.97)	0.98***(0.98, 0.98)	1.00(0.99, 1.01)
Type of care, when time from symptom onset to care=0(Comparison = Outpatient)	21.57***(20.70, 22.48)	43.62***(41.64, 45.70)	67.95***(61.32, 75.30)	73.20***(62.93, 85.15)
Akaike information criterion	2,038,299.00	2,312,173.70	881,563.05	198,461.41

CI, confidence interval

For each of the waves, the adjusted HR for the indigenous population compared to the non-indigenous population was greater than 1 and was statistically significant for the first and third waves, indicating a higher hazard of COVID-19 mortality among indigenous populations compared to the non-indigenous counterparts (Table 5). We found that the unadjusted HR in 2020 was 1.69 (95% CI: 1.62, 1.77) compared to 1.62 (95% CI: 1.54, 1.70) in 2021. Similarly, in the adjusted model (model 6), the HR for 2020 decreased to 1.10 (95% CI: 1.05, 1.15) and the HR for 2021 decreased to 1.10 (95% CI: 1.04, 1.15).

The survival curves for indigenous and non-indigenous populations, stratified by hospitalized and outpatient status in Mexico are shown in Figure 5 . The curves reflect a lower survival probability among hospitalized indigenous populations compared to the hospitalized non-indigenous population.

Figure 5

Survival curves for hospitalized and outpatient COVID-19 cases among indigenous and non-indigenous groups.

Discussion

The COVID-19 pandemic has exacerbated the existing health and socio-economic disparities among vulnerable indigenous populations globally because of the systemic structural disadvantages and higher burden of lockdown measures that put these groups at higher risk of infection and death (Carethers, 2021; Katikireddi ). This study assessed the mortality rate per 1000 person-weeks for indigenous and non-indigenous populations in Mexico from February 2020 to March 2022. At the national level, we estimated the COVID-19 mortality rate among indigenous groups to be substantially higher compared to non-indigenous populations (3.25 vs. 1.94 per 1000 person-weeks). Indigenous individuals experienced significant delays in care seeking for COVID-19 infection and a lower survival probability compared to non-indigenous counterparts.

We found that the COVID-19 mortality rate among indigenous populations was 68% higher than the non-indigenous populations (RR = 1.68). This estimate is similar to that reported in an earlier study from Mexico (RR = 1.65) covering the first 5 months of the pandemic (Argoty-Pantoja ). Another study from Mexico based on data as of November 2020 also reported a higher case fatality ratio of the indigenous group compared to the non-indigenous group (RR = 1.11) (Mallard ). Increased vulnerability and mortality among indigenous populations as compared to the general population during past pandemics is well established (CDC, 2009; Dahal ; Kelm, 1999; Wilson ). For example, during the 1918 influenza pandemic, the death rate among the Māori population was at least 7.3 times higher than that of the population of European descent in New Zealand (Wilson ). Also, during the 2009 H1N1 influenza pandemic, indigenous populations from Australia, Canada, and New Zealand had 3-8 times higher rates of hospitalization and death compared to non-indigenous populations (CDC, 2009).

We observed the highest COVID-19 mortality rate during the first wave followed by the third wave both for indigenous and non-indigenous groups. Indigenous populations were more heavily affected than the non-indigenous consistently across the waves. Social distancing guidelines were not strictly in place and vaccination was not yet available during the first wave whereas the third wave included the period during which the more contagious, severe, and deadly Delta variant of the Coronavirus was circulating (CDC, 2021). RR increased from 1.55 in first wave to 1.65 in second wave to 1.86 in third wave and 2.40 in fourth wave, indicating that over subsequent pandemic waves, the indigenous groups were more likely to succumb to COVID-19. Similarly, the increase in RR from 1.67 in 2020 to 1.75 in 2021 could reflect the effect of social disparities in adopting non-pharmaceutical interventions (NPIs) such as handwashing, using facemasks, less frequent use of public transportation, and social distancing as well as vaccine uptake (Mamelund ). Among the indigenous Mexican populations, precarious socio-economic conditions and factors such as return of the indigenous population to their communities due to the pandemic, poor access to water, language barriers, and limited access to the Internet are important social factors that affect the adoption of preventive measures against COVID-19 (Díaz de León-Martínez ). However, no such data that compares the uptake of NPIs and vaccination among indigenous and non-indigenous groups is available for Mexico. The COVID-19 pandemic started as a rich man's disease within affluent communities and then spread to the poorer section, producing more severe outcomes (Bengali ; Khlat and Le Coeur, 2021; Plümper and Neumayer, 2020). This is consistent with indigenous groups exhibiting worsening mortality outcomes over subsequent pandemic waves. We also found that the R for both indigenous and non-indigenous groups was highest during the fourth wave (1.4 among indigenous vs. 1.5 among non-indigenous). Overall, the R for indigenous and non-indigenous groups was comparable, possibly reflecting the fact that indigenous groups in Mexico are well connected to the rest of the population through internal migration. For example, indigenous population in Mexico migrates to urban destinations where they occupy jobs that increase their vulnerability to COVID-19 infection and they in turn spread the disease to their home communities upon return (Díaz de León-Martínez ). However, the increased RR value could also indicate that the improvements in care and treatment of severe cases became more available to affluent patients than to the indigenous groups as the pandemic progressed.

COVID-19 vaccine uptake among indigenous populations may be influenced by several factors including the availability of accurate, accessible, and culturally relevant information, poor access to the vaccination sites, mistrust in government services, concerns over vaccine safety, past experiences of racism and abuse in healthcare settings, and the lack of involvement of indigenous communities in the planning and implementation of vaccination programs (Castillo ). Many indigenous people live in geographically rural or remote areas, inhabit multigenerational housing, are poor, face food insecurity, have increased rates of chronic diseases, and have limited access to running water. Consequently, the one-size-fits-all response that was frequently implemented to address COVID-19 was often ineffective for these vulnerable groups (Power ; United Nations Department of Economic and Social Affairs, 2020). In Mexico specifically, indigenous populations have a higher prevalence of metabolic syndrome and its components such as low HDL-cholesterol levels, central obesity and elevated blood pressure (Mendoza-Caamal ) that increase the risk of severe COVID-19 (Steenblock ). In addition, indigenous communities in Mexico suffer from higher levels of social deprivation including lagging in education, access to health services and social security, housing quality and space, and basic housing services (Díaz de León-Martínez ; Consejo Nacional de Evaluación de la Politica de Desrrollo Social, 2019). As in many countries, the Mexican health administration recommended COVID-19 vaccination following sequential prioritization to health workers, people aged 60 years and above, and people aged 50-59 with comorbidities followed by the general population (de Vacunación Covid, 2021). However, social and ethnic vulnerability including indigenous status was not taken into account when designing these recommendations. Indeed, a 2017 review of pandemic preparedness plans across nations revealed that the risk groups for vaccination were based on medical conditions without regard to social risks (Mamelund, 2017; Mamelund and Dimka, 2021).

In Mexico, 23 federal entities out of 32 had a higher mortality rate among indigenous compared to non-indigenous, and this difference was statistically significant in 13 entities. For 9 states including Chiapas, Durango, and Nuevo Leon, the mortality rate was higher among non-indigenous groups, though not statistically significant. Further studies are needed to better understand the drivers leading to higher mortality among non-indigenous populations in these areas.

In summary, we found indigenous status to be an important risk factor for COVID-19 mortality in Mexico. While the crude hazard of death among the indigenous population was 67.1% higher than that of the non-indigenous, the hazard of death was only 8.0% greater among indigenous when controlling for sex, age, comorbidity category, and the interaction between type of care and care-seeking delay. This indicates that ‘indigenous status’ alone is not a significant factor contributing to a higher risk of deaths due to COVID-19 among the indigenous population, but the associated socio-demographic, cultural and care-seeking practices in this group explain such differences. Moreover, epigenetic mechanisms could also mediate the effect of racialized experiences of colonization, oppression, and violence on COVID-19 outcomes among indigenous populations (Curtice and Choo, 2020; Encinosa, 2021). Similar to our Mexican study, in Brazil, which has high ethnic diversity, Pardo ethnicity was found to be the second most important risk factor for death due to COVID-19 after age. Compared to White Brazilians, Pardo and Black Brazilians who were admitted to the hospital suffered a significantly higher risk of mortality (HR for Pardo: 1.45; 95% CI: 1.33, 1.58) (Baqui ). Consistent with this study, we found declining survival curves for Mexican indigenous groups compared to non-indigenous for both the outpatient and hospitalized COVID-19 positive cases.

We found a significantly higher proportion of deaths among total cases in the indigenous population based on hospitalization and outpatient status. This was consistent with previous studies from Mexico (Argoty-Pantoja ; Ibarra-Nava ). We also found a statistically significant higher average care-seeking delay for the indigenous population compared to non-indigenous (4.41 days vs. 4.07 days) which could be due in part to longer travel distances between place of residence and health facility for the indigenous population compared to the non-indigenous populations besides other socio-economic disparities affecting health service utilization.

Our study has some notable limitations. First, the nature of available variables in the dataset did not allow us to fully understand the drivers behind the higher mortality rates in indigenous populations compared to the general population. For instance, variables related to COVID-19 vaccination status, self-reported adherence to NPIs such as mask wearing and shelter-in-place orders, and access to health care services could help us better understand COVID-19 mortality disparities. Results for waves other than the first wave may be confounded by the disparity in vaccine uptake and disparity in NPI use. It is also worth noting the possibility of underreporting of COVID-19 deaths among indigenous populations as reported for Brazil (Fellows ) where it was estimated a 103% underreporting of COVID-19 deaths among indigenous populations. Similarly, according to a previous study from Mexico (Dahal ), confirmed COVID-19 deaths accounted for only ∼38% of total excess deaths estimated for Mexico partly due to underreporting. Our results are therefore underestimating the toll of the disease burden on the indigenous peoples. Differences in underreporting between indigenous and non-indigenous populations need further investigation.

Conclusion

Our study found that the COVID-19 mortality rate among indigenous populations was 68% higher than that for the non-indigenous population (RR = 1.68) from February 2020 to March 2022, and both the R and RR were highest during the fourth wave of the COVID-19 pandemic, reflecting the possible impact of vaccine uptake and NPI use disparities between indigenous and non-indigenous population. Our findings indicate that indigenous status is an important risk factor for mortality in Mexico such that the hazard of death among indigenous populations was 67% higher among indigenous people compared to non-indigenous groups in the unadjusted model and 8% greater when controlling for sex, age, number of comorbidities and the interaction between type of care and care-seeking delay. Indigenous people showed higher delays in care seeking for COVID-19 infection and a lower survival probability compared to non-indigenous populations, pointing to a link between care seeking practices and indigenous status. Further studies are warranted to disentangle the mechanisms through which indigenous populations are disproportionately affected by COVID-19 compared to non-indigenous populations.

Author contributions

GC and SD conceptualized the study; SD, SEM, RL, LS, SSB, and GC designed the methodology; SD prepared the original draft; SD performed the formal analysis; and SEM, RL, LS, SSB, GC edited and reviewed the manuscript., All the authors read and approved the final version of this manuscript.

Data availability

All the data used in the study are publicly available. The sources of those data are provided under the methods section.

Funding

SD was funded by 2CI Doctoral Fellowship at Georgia State University. This research did not receive any specific grant from funding agencies in the public, commercial, or not-for profit sectors.

Declaration of Competing Interest

The authors have no competing interests to declare.