How can I keep the order of the row of dataframe while cbinding?

I am trying to make a function for making a neat table from the result of HSD.test from the agricolae package. The output of the HSD.test is a list where the treatment lettering is sorted in descending order of the mean. I would like to arrange "mean ± sd letter' by T1, T2,... so on. Though mean and sd are arranged this way, the lettering does not follow the order. I cannot accomplish the following with the following function, as the letters are not aligned properly with the treatment (which is the row names in each data frame below). I have searched here and there and tried to collected code fragments to make the function so far.

How can I keep the order of the row of df3 during cbinding?

MyDf<-data.frame(treatment = rep(c('T1','T2','T3','T4'), each = 3),
                 p1 = c(28.5, 21.7, 23.0, 14.9, 10.6, 13.1, 41.8, 39.2, 28.0, 38.2, 40.4, 32.1), 
                 p2 = c(32, 37, 36, 23, 28, 22, 67, 52, 55, 18, 27, 17),
                 p3 = c(5.6, 3.7, 4.9, 7.1, 11.3, 14, 2.3, 5.4, 3.3, 11.6, 10.1, 12)
                 )

MeltMyDf<-melt(MyDf)

MyDfSE<-summarySE(MeltMyDf, measurevar = "value", groupvars = c("variable", "treatment"))

x <- as.character(unique(MyDfSE$variable))

MyModels<-sapply(x, function(my) {lm(value~treatment, data=MeltMyDf, variable==my)}, simplify=FALSE)

MyGroups<- lapply(MyModels, function(m) HSD.test((m), "treatment", alpha = 0.05, group = TRUE, console = FALSE, variable==variable))

#---- the function I have created ------------
make_HSD_table<-function(HSDlist){
  
  mycolnames<-names(HSDlist)
  mycolnumber<-length(mycolnames)
  
  df1<-as.data.frame(lapply(lapply(HSDlist, `[[`, 'means'), '[', 'value'))
  df1[order(row.names(df1)), ]
  df1<-round(df1, digits = 2)
  
  df2<-as.data.frame(lapply(lapply(HSDlist, `[[`, 'means'), '[', c('std')))
  df2[order(row.names(df2)), ]
  df2<-round(df2, digits = 2)
  
  df3<-as.data.frame(lapply(lapply(HSDlist, `[[`, 'groups'), '[', 'groups'))
  df3[order(row.names(df3)), ]

  myrownumber<-nrow(df1)
  pm_df<-data.frame(replicate(mycolnumber, rep('±', myrownumber)))
  
  my_table <- cbind(df1, pm_df, df2, df3)[order(c(seq_along(df1), seq_along(pm_df),seq_along(df2), seq_along(df3)))]
  
  
  my_table <- cbind(sapply(split.default(my_table, as.integer(gl(ncol(my_table), 4, ncol(my_table)))), function(x) do.call(paste, x)))
  
  colnames(my_table)<-mycolnames
  rownames(my_table)<-rownames(df1)

  write.table(my_table, 'my_HSD_table.csv', append = F, sep = ',', row.names = FALSE)
  
  return(my_table)
}

#-----------------------------------------
make_HSD_table(MyGroups)

   p1                p2               p3             
T1 "24.4 ± 3.61 a"   "35 ± 2.65 a"    "4.73 ± 0.96 a"
T2 "12.87 ± 2.16 ab" "24.33 ± 3.21 b" "10.8 ± 3.48 a"
T3 "36.33 ± 7.33 bc" "58 ± 7.94 bc"   "3.67 ± 1.58 b"
T4 "36.9 ± 4.3 c"    "20.67 ± 5.51 c" "11.23 ± 1 b"  

**As you can see, the values±sd are sorted according to the treatment. But the letters are not sorted and placed with wrong mean±sd!**

Finally, I could correct my function to make an output of mean ± sd grouping_letter table from HSD.test which is sorted in the order of data frame. Here is the workflow:

library(reshape2)
library(Rmisc)
library(agricolae)
library(dplyr)
library(stringr)
MyDf<-data.frame(treatment = rep(c('T1','T2','T3','T4'), each = 3),
p1 = c(28.5, 21.7, 23.0, 14.9, 10.6, 13.1, 41.8, 39.2, 28.0, 38.2, 40.4, 32.1), 
p2 = c(32, 37, 36, 23, 28, 22, 67, 52, 55, 18, 27, 17),
p3 = c(5.6, 3.7, 4.9, 7.1, 11.3, 14, 2.3, 5.4, 3.3, 11.6, 10.1, 12)
)
MeltMyDf<-melt(MyDf)
MyDfSE<-summarySE(MeltMyDf, measurevar = "value", groupvars = c("variable", "treatment"))
x <- as.character(unique(MyDfSE$variable))
MyModels<-sapply(x, function(my) {lm(value~treatment, data=MeltMyDf, variable==my)}, simplify=FALSE)
MyGroups<- lapply(MyModels, function(m) HSD.test((m), "treatment", alpha = 0.05, group = TRUE, console = FALSE, variable==variable))

Now, the function for calling the table. It will also export the table as a csv file.

make_HSD_table<-function(HSDlist){
mycolnames<-names(HSDlist)
mycolnumber<-length(mycolnames)
df1<-as.data.frame(lapply(lapply(HSDlist, `[[`, 'means'), '[', 'value'))
df1[order(row.names(df1)), ]
df1<-round(df1, digits = 2)
df2<-as.data.frame(lapply(lapply(HSDlist, `[[`, 'means'), '[', c('std')))
df2[order(row.names(df2)), ]
df2<-round(df2, digits = 2)
df3<-lapply(lapply(HSDlist, `[[`, 'groups'), '[', 'groups')
df3<-Map(cbind, df3, my_row_names = lapply(df3, rownames))
df3<-lapply(df3, function(df) {df[order(df$my_row_names), ]})
df3<-lapply(df3, function(x) x[,1])
myrownumber<-nrow(df1)
pm_df<-data.frame(replicate(mycolnumber, rep('±', myrownumber)))
my_table <- cbind(df1, pm_df, df2, df3)[order(c(seq_along(df1), seq_along(pm_df),seq_along(df2), seq_along(df3)))]
my_table <- cbind(sapply(split.default(my_table, as.integer(gl(ncol(my_table), 4, ncol(my_table)))), function(x) do.call(paste, x)))
colnames(my_table)<-mycolnames
my_table<-data.frame(treatment = rownames(df1), my_table)
write.table(my_table, 'my_HSD_table.csv', append = F, sep = ',', row.names = FALSE)
return(my_table)
}

Let's make the table from HSD.test output.

make_HSD_table(MyGroups)
treatment              p1              p2            p3
1        T1  24.4 ± 3.61 bc     35 ± 2.65 b 4.73 ± 0.96 b
2        T2  12.87 ± 2.16 c 24.33 ± 3.21 bc 10.8 ± 3.48 a
3        T3 36.33 ± 7.33 ab     58 ± 7.94 a 3.67 ± 1.58 b
4        T4    36.9 ± 4.3 a  20.67 ± 5.51 c   11.23 ± 1 a

As you can see, the letters have been properly assigned to the respective mean±sd columns.