I have performed one way ANOVA's on different categories of my data using R. But I am confused how to perform Tukey's test for each category. My data looks like:

my_data<-
   Entry       gene        expression
     Abx    A0A2N8I1K7       1.11
    Abx    A0A2N8I1K7       1.11
    Abx    A0A2N8I1K7       1.11
    Afx     A0A2N8I1K7       1.10
    Afx    A0A2N8I1K7       1.10
    Afx    A0A2N8I1K7       1.10
     B     A0A2N8I1K7       1.08
      B    A0A2N8I1K7       1.08
      B     A0A2N8I1K7       1.08
      B     A0A2N8I1K7       1.08
    BY     A0A2N8I1K7       1.09
     BY    A0A2N8I1K7       1.10
    BY    A0A2N8I1K7       1.08
    Abx    A0A2N8I321       1.03
    Abx   A0A2N8I321       1.03
    Abx   A0A2N8I321       1.03
     Afx   A0A2N8I321       1.01
     Afx   A0A2N8I321       1.01
     Afx   A0A2N8I321       1.01
      B     A0A2N8I321       0.98
      B    A0A2N8I321       1.00
      B    A0A2N8I321       0.97
     B   A0A2N8I321       0.98
     BY  A0A2N8I321       1.02
     BY  A0A2N8I321       1.01
     BY  A0A2N8I321       1.00
    Abx  A0A2N8I3R9       1.02
    Abx A0A2N8I3R9       1.02
    Abx  A0A2N8I3R9       1.01
    Afx A0A2N8I3R9       1.00
   Afx A0A2N8I3R9       1.01
32    Afx A0A2N8I3R9       1.00
33      B A0A2N8I3R9       0.98
34      B A0A2N8I3R9       0.98
35      B A0A2N8I3R9       0.98
36      B A0A2N8I3R9       0.98
37     BY A0A2N8I3R9       0.99
38     BY A0A2N8I3R9       0.99
39     BY A0A2N8I3R9       0.99

There are more genes in my data. I have performed one-way ANOVA for each gene or category using following formula:

a<-lapply(split(my_data, my_data$gene), function(i){
  aov(lm(expression ~ Entry, data = i))
})

Then I did Tukey's test for example, for gene, A0A2N8I1K7 like this:

TukeyHSD(a$A0A2N8I1K7)

But this procedure for doing Tukey's test is so time-consuming as I have lots of genes or categories in my data. Is there any other way to perform Tukey's test?